diff --git a/Detector/DetCRD/CMakeLists.txt b/Detector/DetCRD/CMakeLists.txt
index 6062572c26a34c19e31d081837ea89cddedd8dd9..177a4672ffa781050400cf8aff82f8a593d13beb 100644
--- a/Detector/DetCRD/CMakeLists.txt
+++ b/Detector/DetCRD/CMakeLists.txt
@@ -22,6 +22,7 @@ gaudi_add_module(DetCRD
src/Tracker/SiTrackerStitching_v01_geo.cpp
src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
src/Tracker/SiTrackerStaggeredLadder_v02_geo.cpp
+ src/Tracker/SiTrackerStaggeredLadder_v03_geo.cpp
src/Tracker/SiTrackerComposite_v01_geo.cpp
src/Tracker/SiTrackerComposite_v02_geo.cpp
src/Tracker/TPC_Simple_o1_v01.cpp
@@ -30,6 +31,8 @@ gaudi_add_module(DetCRD
src/Tracker/SiTracker_itkbarrel_v02_geo.cpp
src/Tracker/SiTracker_otkbarrel_v01_geo.cpp
src/Tracker/SiTracker_otkendcap_v01_geo.cpp
+ src/Tracker/SiTracker_otkbarrel_v02_geo.cpp
+ src/Tracker/SiTracker_otkendcap_v02_geo.cpp
src/Other/ParaffinEndcap_v01.cpp
LINK ${DD4hep_COMPONENT_LIBRARIES}
diff --git a/Detector/DetCRD/compact/CRD_common_v02/ITK_StaggeredStave_v03_01.xml b/Detector/DetCRD/compact/CRD_common_v02/ITK_StaggeredStave_v03_01.xml
new file mode 100644
index 0000000000000000000000000000000000000000..c6bbe28055d078eb7429563b5cf519e18136a66c
--- /dev/null
+++ b/Detector/DetCRD/compact/CRD_common_v02/ITK_StaggeredStave_v03_01.xml
@@ -0,0 +1,126 @@
+<lccdd>
+ <define>
+ <constant name="ITKBarrel_sensitive_thickness" value="0.15*mm"/>
+ <constant name="ITKBarrel_inner_radius_1" value="ITKBarrel1_inner_radius"/>
+ <constant name="ITKBarrel_inner_radius_2" value="ITKBarrel2_inner_radius"/>
+ <constant name="ITKBarrel_inner_radius_3" value="ITKBarrel3_inner_radius"/>
+ <constant name="ITKBarrel_ladder_length_1" value="ITKBarrel1_half_length*2"/>
+ <constant name="ITKBarrel_ladder_length_2" value="ITKBarrel2_half_length*2"/>
+ <constant name="ITKBarrel_ladder_length_3" value="ITKBarrel3_half_length*2"/>
+ <constant name="ITKBarrel_sensor_size" value="20*mm"/>
+ <constant name="ITKBarrel_sensor_side" value="0.3*mm"/>
+ <constant name="ITKBarrel_sensor_dead" value="2.0*mm"/>
+ <constant name="ITKBarrel_sensor_gap" value="0.1*mm"/>
+ <constant name="ITKBarrel_sensor_number" value="14"/>
+ <constant name="ITKBarrel_module_length" value="(ITKBarrel_sensor_size+ITKBarrel_sensor_gap)*ITKBarrel_sensor_number/2
+ -ITKBarrel_sensor_gap"/>
+ <constant name="ITKBarrel_module_width" value="2*ITKBarrel_sensor_size+ITKBarrel_sensor_gap"/>
+ <constant name="ITKBarrel_module_gap" value="0.3*mm"/>
+ <constant name="ITKBarrel_module_number_1" value="7"/>
+ <constant name="ITKBarrel_module_number_2" value="10"/>
+ <constant name="ITKBarrel_module_number_3" value="14"/>
+ </define>
+
+ <detectors>
+ <detector id="DetID_ITKBarrel" name="ITKBarrel" type="SiTrackerStaggeredLadder_v03" vis="SeeThrough"
+ readout="ITKBarrelCollection" combineHits="true" printLevel="INFO" insideTrackingVolume="true">
+ <envelope>
+ <shape type="Assembly"/>
+ </envelope>
+
+ <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_PIXEL "/>
+
+ <layer layer_id="0" phi0="0" n_ladders="44" radius="ITKBarrel_inner_radius_1" rotate="-55.85*mrad">
+ <support length="ITKBarrel_ladder_length_1" width="ITKBarrel_module_width" material="Air" vis="SeeThrough">
+ <slice name="TrussFrame" thickness="208*um" width="ITKBarrel_module_width" material="CF_ITK" vis="LightGrayVis"/>
+ <slice name="CarbonFleece" thickness=" 20*um" width="ITKBarrel_module_width" material="CarbonFleece_ITK" vis="LightGrayVis"/>
+ <slice name="GraphiteFoil" thickness=" 30*um" width="ITKBarrel_module_width" material="Graphite_ITK" vis="GrayVis"/>
+ <slice name="CoolingTube" thickness=" 64*um" width="ITKBarrel_module_width" material="Polyimide_ITK" vis="SeeThrough"/>
+ <slice name="CoolingFluid" thickness="190*um" width="ITKBarrel_module_width" material="EquivalWater_ITK" vis="SeeThrough"/>
+ <slice name="CFPlate" thickness="150*um" width="ITKBarrel_module_width" material="CF_ITK" vis="GrayVis"/>
+ <slice name="CarbonFleece" thickness=" 20*um" width="ITKBarrel_module_width" material="CarbonFleece_ITK" vis="LightGrayVis"/>
+ <slice name="Glue" thickness="100*um" width="ITKBarrel_module_width" material="CER_ITK" vis="GrayVis"/>
+ </support>
+ <sensitive length="ITKBarrel_ladder_length_1" width="ITKBarrel_module_width" thickness="ITKBarrel_sensitive_thickness" material="G4_Si" vis="OrangeVis">
+ <module length="ITKBarrel_module_length" width="ITKBarrel_module_width" material="G4_Si" vis="CyanVis">
+ <row repeat="1" gap="0" up_side="0"/>
+ <column repeat="ITKBarrel_module_number_1" gap="ITKBarrel_module_gap" left_side="0"/>
+ <sensor repeat="ITKBarrel_sensor_number" length="ITKBarrel_sensor_size-2*ITKBarrel_sensor_side" width="ITKBarrel_sensor_size-2*ITKBarrel_sensor_side-ITKBarrel_sensor_dead"
+ gap="ITKBarrel_sensor_gap" dead="ITKBarrel_sensor_dead" side="ITKBarrel_sensor_side" material="G4_Si" vis="MagentaVis"/>
+ </module>
+ </sensitive>
+ <flex length="ITKBarrel_ladder_length_1" width="ITKBarrel_module_width" material="Air" vis="SeeThrough">
+ <slice name="Glue" thickness="100*um" material="CER_ITK" vis="YellowVis"/>
+ <slice name="FPCInsulating" thickness="100*um" material="Polyimide_ITK" vis="YellowVis"/>
+ <slice name="FPCMetal" thickness="100*um" material="G4_Al" vis="GrayVis"/>
+ <slice name="OEComponent1" thickness=" 25*um" material="Kapton" vis="YellowVis"/>
+ <slice name="OEComponent2" thickness=" 56*um" material="G4_POLYETHYLENE" vis="GreenVis"/>
+ <slice name="OEComponent3" thickness=" 3*um" material="G4_Cu" vis="RedVis"/>
+ </flex>
+ </layer>
+
+ <layer layer_id="1" phi0="0" n_ladders="64" radius="ITKBarrel_inner_radius_2" rotate="-55.85*mrad">
+ <support length="ITKBarrel_ladder_length_2" width="ITKBarrel_module_width" material="Air" vis="SeeThrough">
+ <slice name="TrussFrame" thickness="208*um" width="ITKBarrel_module_width" material="CF_ITK" vis="LightGrayVis"/>
+ <slice name="CarbonFleece" thickness=" 20*um" width="ITKBarrel_module_width" material="CarbonFleece_ITK" vis="LightGrayVis"/>
+ <slice name="GraphiteFoil" thickness=" 30*um" width="ITKBarrel_module_width" material="Graphite_ITK" vis="GrayVis"/>
+ <slice name="CoolingTube" thickness=" 64*um" width="ITKBarrel_module_width" material="Polyimide_ITK" vis="SeeThrough"/>
+ <slice name="CoolingFluid" thickness="190*um" width="ITKBarrel_module_width" material="EquivalWater_ITK" vis="SeeThrough"/>
+ <slice name="CFPlate" thickness="150*um" width="ITKBarrel_module_width" material="CF_ITK" vis="GrayVis"/>
+ <slice name="CarbonFleece" thickness=" 20*um" width="ITKBarrel_module_width" material="CarbonFleece_ITK" vis="LightGrayVis"/>
+ <slice name="Glue" thickness="100*um" width="ITKBarrel_module_width" material="CER_ITK" vis="GrayVis"/>
+ </support>
+ <sensitive length="ITKBarrel_ladder_length_2" width="ITKBarrel_module_width" thickness="ITKBarrel_sensitive_thickness" gap="0.3*mm" material="G4_Si" vis="OrangeVis">
+ <module length="ITKBarrel_module_length" width="ITKBarrel_module_width" material="G4_Si" vis="CyanVis">
+ <row repeat="1" gap="0" up_side="0"/>
+ <column repeat="ITKBarrel_module_number_2" gap="ITKBarrel_module_gap" left_side="0"/>
+ <sensor repeat="ITKBarrel_sensor_number" length="ITKBarrel_sensor_size-2*ITKBarrel_sensor_side" width="ITKBarrel_sensor_size-ITKBarrel_sensor_dead"
+ gap="ITKBarrel_sensor_gap" dead="ITKBarrel_sensor_dead" side="ITKBarrel_sensor_side" material="G4_Si" vis="MagentaVis"/>
+ </module>
+ </sensitive>
+ <flex length="ITKBarrel_ladder_length_2" width="ITKBarrel_module_width" material="Air" vis="SeeThrough">
+ <slice name="Glue" thickness="100*um" material="CER_ITK" vis="YellowVis"/>
+ <slice name="FPCInsulating" thickness="100*um" material="Polyimide_ITK" vis="YellowVis"/>
+ <slice name="FPCMetal" thickness="100*um" material="G4_Al" vis="GrayVis"/>
+ <slice name="OEComponent1" thickness=" 25*um" material="Kapton" vis="YellowVis"/>
+ <slice name="OEComponent2" thickness=" 56*um" material="G4_POLYETHYLENE" vis="GreenVis"/>
+ <slice name="OEComponent3" thickness=" 3*um" material="G4_Cu" vis="RedVis"/>
+ </flex>
+ </layer>
+
+ <layer layer_id="2" phi0="0" n_ladders="102" radius="ITKBarrel_inner_radius_3" rotate="-55.85*mrad">
+ <support length="ITKBarrel_ladder_length_3" width="ITKBarrel_module_width" material="Air" vis="SeeThrough">
+ <slice name="TrussFrame" thickness="208*um" width="ITKBarrel_module_width" material="CF_ITK" vis="LightGrayVis"/>
+ <slice name="CarbonFleece" thickness=" 20*um" width="ITKBarrel_module_width" material="CarbonFleece_ITK" vis="LightGrayVis"/>
+ <slice name="GraphiteFoil" thickness=" 30*um" width="ITKBarrel_module_width" material="Graphite_ITK" vis="GrayVis"/>
+ <slice name="CoolingTube" thickness=" 64*um" width="ITKBarrel_module_width" material="Polyimide_ITK" vis="SeeThrough"/>
+ <slice name="CoolingFluid" thickness="190*um" width="ITKBarrel_module_width" material="EquivalWater_ITK" vis="SeeThrough"/>
+ <slice name="CFPlate" thickness="150*um" width="ITKBarrel_module_width" material="CF_ITK" vis="GrayVis"/>
+ <slice name="CarbonFleece" thickness=" 20*um" width="ITKBarrel_module_width" material="CarbonFleece_ITK" vis="LightGrayVis"/>
+ <slice name="Glue" thickness="100*um" width="ITKBarrel_module_width" material="CER_ITK" vis="GrayVis"/>
+ </support>
+ <sensitive length="ITKBarrel_ladder_length_3" width="ITKBarrel_module_width" thickness="ITKBarrel_sensitive_thickness" gap="0.3*mm" material="G4_Si" vis="OrangeVis">
+ <module length="ITKBarrel_module_length" width="ITKBarrel_module_width" material="G4_Si" vis="CyanVis">
+ <row repeat="1" gap="0" up_side="0"/>
+ <column repeat="ITKBarrel_module_number_3" gap="ITKBarrel_module_gap" left_side="0"/>
+ <sensor repeat="ITKBarrel_sensor_number" length="ITKBarrel_sensor_size-2*ITKBarrel_sensor_side" width="ITKBarrel_sensor_size-ITKBarrel_sensor_dead"
+ gap="ITKBarrel_sensor_gap" dead="ITKBarrel_sensor_dead" side="ITKBarrel_sensor_side" material="G4_Si" vis="MagentaVis"/>
+ </module>
+ </sensitive>
+ <flex length="ITKBarrel_ladder_length_3" width="ITKBarrel_module_width" material="Air" vis="SeeThrough">
+ <slice name="Glue" thickness="100*um" material="CER_ITK" vis="YellowVis"/>
+ <slice name="FPCInsulating" thickness="100*um" material="Polyimide_ITK" vis="YellowVis"/>
+ <slice name="FPCMetal" thickness="100*um" material="G4_Al" vis="GrayVis"/>
+ <slice name="OEComponent1" thickness=" 25*um" material="Kapton" vis="YellowVis"/>
+ <slice name="OEComponent2" thickness=" 56*um" material="G4_POLYETHYLENE" vis="GreenVis"/>
+ <slice name="OEComponent3" thickness=" 3*um" material="G4_Cu" vis="RedVis"/>
+ </flex>
+ </layer>
+ </detector>
+ </detectors>
+ <readouts>
+ <readout name="ITKBarrelCollection">
+ <id>system:5,side:-2,layer:9,stave:8,module:8,sensor:5</id>
+ </readout>
+ </readouts>
+</lccdd>
diff --git a/Detector/DetCRD/compact/CRD_common_v02/SIT_StaggeredStave_v02.xml b/Detector/DetCRD/compact/CRD_common_v02/SIT_StaggeredStave_v02.xml
index 4329abb09aa5800daace93ac185206d5a6200356..8f8d353d7a67b0ff2efc717fb87ce96860d08e38 100644
--- a/Detector/DetCRD/compact/CRD_common_v02/SIT_StaggeredStave_v02.xml
+++ b/Detector/DetCRD/compact/CRD_common_v02/SIT_StaggeredStave_v02.xml
@@ -4,12 +4,12 @@
<constant name="SIT_sensitive_thickness" value="0.15*mm"/>
<!--constant name="SIT_support_thickness" value="4*mm"/--> <!--use the barestave num-->
<constant name="SIT_module_length" value="140*mm"/> <!--module length-->
- <constant name="SIT_inner_radius_1" value="SIT1_inner_radius"/>
- <constant name="SIT_inner_radius_2" value="SIT2_inner_radius"/>
- <constant name="SIT_inner_radius_3" value="SIT3_inner_radius"/>
- <constant name="SIT_half_length_1" value="SIT1_half_length"/>
- <constant name="SIT_half_length_2" value="SIT2_half_length"/>
- <constant name="SIT_half_length_3" value="SIT3_half_length"/>
+ <constant name="SIT_inner_radius_1" value="ITKBarrel1_inner_radius"/>
+ <constant name="SIT_inner_radius_2" value="ITKBarrel2_inner_radius"/>
+ <constant name="SIT_inner_radius_3" value="ITKBarrel3_inner_radius"/>
+ <constant name="SIT_half_length_1" value="ITKBarrel1_half_length"/>
+ <constant name="SIT_half_length_2" value="ITKBarrel2_half_length"/>
+ <constant name="SIT_half_length_3" value="ITKBarrel3_half_length"/>
</define>
<detectors>
diff --git a/Detector/DetCRD/compact/CRD_common_v02/materials.xml b/Detector/DetCRD/compact/CRD_common_v02/materials.xml
index 868157a2ffc7eef086c5e67b8ea27dd7630717ff..e30604a0373f7e3e8f6cfbd691807c960e5e8297 100644
--- a/Detector/DetCRD/compact/CRD_common_v02/materials.xml
+++ b/Detector/DetCRD/compact/CRD_common_v02/materials.xml
@@ -827,12 +827,46 @@
<fraction n="0.895" ref="Al"/>
</material>
- <material name="ParaffinWax">
- <D value="0.9" unit="g/cm3"/>
- <composite n="25" ref="C"/>
- <composite n="52" ref="H"/>
- </material>
+ <material name="ParaffinWax">
+ <D value="0.9" unit="g/cm3"/>
+ <composite n="25" ref="C"/>
+ <composite n="52" ref="H"/>
+ </material>
+
+ <!--ITK-->
+ <material name="CER_ITK">
+ <D type="density" value="0.880906" unit="g/cm3"/>
+ <composite n="6" ref="C"/>
+ <composite n="5" ref="H"/>
+ <composite n="3" ref="N"/>
+ <composite n="3" ref="O"/>
+ </material>
+
+ <material name="CarbonFleece_ITK">
+ <D type="density" value="0.3997838" unit="g/cm3" />
+ <composite n="1" ref="C" />
+ </material>
+
+ <material name="Graphite_ITK">
+ <D type="density" value="1.607564" unit="g/cm3" />
+ <composite n="1" ref="C" />
+ </material>
+ <material name="CF_ITK">
+ <D type="density" value="1.60885" unit="g/cm3"/>
+ <fraction n="1.0" ref="CarbonFiber"/>
+ </material>
+
+ <material name="Polyimide_ITK">
+ <D type="density" value="1.428233" unit="g/cm3"/>
+ <fraction n="1.0" ref="Polyimide"/>
+ </material>
+
+ <!-- to reduce equivalent thickness -->
+ <material name="EquivalWater_ITK">
+ <D type="density" value="2" unit="g/cm3"/>
+ <fraction n="1.0" ref="G4_WATER"/>
+ </material>
</materials>
<surfaces>
diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml
index ea9d8d90df5fd06e5423fa9809ae8094cd03ade0..2f0a6965581f8694c93947c66f011d41314230e0 100644
--- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml
+++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml
@@ -113,12 +113,12 @@
<constant name="TPC_half_length" value="2900*mm"/>
<constant name="OuterTracker_half_length" value="TPC_half_length"/>
- <constant name="SIT1_inner_radius" value="240*mm"/>
- <constant name="SIT2_inner_radius" value="350*mm"/>
- <constant name="SIT3_inner_radius" value="(600-30)*mm"/>
- <constant name="SIT1_half_length" value="500.5*mm"/>
- <constant name="SIT2_half_length" value="715*mm"/>
- <constant name="SIT3_half_length" value="1001*mm"/>
+ <constant name="ITKBarrel1_inner_radius" value="235*mm"/>
+ <constant name="ITKBarrel2_inner_radius" value="345*mm"/>
+ <constant name="ITKBarrel3_inner_radius" value="555.6*mm"/>
+ <constant name="ITKBarrel1_half_length" value="500.5*mm"/>
+ <constant name="ITKBarrel2_half_length" value="715*mm"/>
+ <constant name="ITKBarrel3_half_length" value="1001*mm"/>
<!-- Parameters of time of flight tracker -->
<constant name="OTKBarrel_inner_radius" value="1800*mm"/>
@@ -135,16 +135,16 @@
<constant name="SiTracker_endcap_barrel_zgap" value="5*mm"/>
<constant name="SiTracker_endcap_barrel_rgap" value="-5*mm"/>
<constant name="SiTracker_endcap_gas_zgap" value="3*mm"/>
- <constant name="SiTracker_endcap_gas_rgap" value="20*mm"/>
- <constant name="SiTracker_endcap_z1" value="SIT1_half_length+SiTracker_endcap_barrel_zgap"/>
- <constant name="SiTracker_endcap_z2" value="SIT2_half_length+SiTracker_endcap_barrel_zgap"/>
- <constant name="SiTracker_endcap_z3" value="SIT3_half_length+SiTracker_endcap_barrel_zgap"/>
+ <constant name="SiTracker_endcap_gas_rgap" value="10*mm"/>
+ <constant name="SiTracker_endcap_z1" value="ITKBarrel1_half_length+SiTracker_endcap_barrel_zgap"/>
+ <constant name="SiTracker_endcap_z2" value="ITKBarrel2_half_length+SiTracker_endcap_barrel_zgap"/>
+ <constant name="SiTracker_endcap_z3" value="ITKBarrel3_half_length+SiTracker_endcap_barrel_zgap"/>
<constant name="SiTracker_endcap_z4" value="1500*mm"/>
<constant name="SiTracker_endcap_z5" value="TPC_half_length+SiTracker_endcap_gas_zgap"/>
- <constant name="SiTracker_endcap_outer_radius1" value="SIT1_inner_radius"/>
- <constant name="SiTracker_endcap_outer_radius2" value="SIT2_inner_radius"/>
- <constant name="SiTracker_endcap_outer_radius3" value="SIT3_inner_radius-SiTracker_endcap_gas_rgap"/>
- <constant name="SiTracker_endcap_outer_radius4" value="SIT3_inner_radius-SiTracker_endcap_gas_rgap"/>
+ <constant name="SiTracker_endcap_outer_radius1" value="ITKBarrel1_inner_radius+SiTracker_endcap_gas_rgap"/>
+ <constant name="SiTracker_endcap_outer_radius2" value="ITKBarrel2_inner_radius+SiTracker_endcap_gas_rgap"/>
+ <constant name="SiTracker_endcap_outer_radius3" value="ITKBarrel3_inner_radius"/>
+ <constant name="SiTracker_endcap_outer_radius4" value="ITKBarrel3_inner_radius"/>
<constant name="SiTracker_endcap_outer_radius5" value="TPC_outer_radius+SiTracker_endcap_barrel_rgap"/>
<constant name="SiTracker_endcap_inner_radius1" value="75*mm"/>
<constant name="SiTracker_endcap_inner_radius2" value="101.9*mm"/>
diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml
index 34a822bccdfc88cd35e1fe55de866efea5477659..85bfb00dc61a5135d49efdc1b43fe8fd10231e24 100644
--- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml
+++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml
@@ -32,7 +32,8 @@
<include ref="../CRD_common_v02/VXD_Composite_v01_02.xml"/>
<include ref="../CRD_common_v02/FTD_SkewRing_v01_07.xml"/>
<!--include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/-->
- <include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/>
+ <!--include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/-->
+ <include ref="../CRD_common_v02/ITK_StaggeredStave_v03_01.xml"/>
<!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/-->
<!-- use 10 rows clustering version-->
<include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/>
diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml
index 4ddc0ae2f2ff11ab9d78900490f4bcf3f8fdee3c..eb449954dc9c77467d22c6659383df79c5811e4c 100644
--- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml
+++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml
@@ -33,7 +33,8 @@
<include ref="../CRD_common_v02/VXD_Composite_v01_02.xml"/>
<include ref="../CRD_common_v02/FTD_SkewRing_v01_07.xml"/>
<!--include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/-->
- <include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/>
+ <!--include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/-->
+ <include ref="../CRD_common_v02/ITK_StaggeredStave_v03_01.xml"/>
<!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/-->
<!--use 10 rows clustering version/-->
<include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/>
diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index 0b9304d9202bc1c54638bfc3d1d5a2ffd843d7fc..6bf8a8e7cf1e59bea6b6bb6abf0ef9f3900425b4 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
@@ -49,9 +49,8 @@ podioinput = PodioInput("PodioReader", collections=[
# "EventHeader",
"MCParticle",
"VXDCollection",
- "SITCollection",
+ "ITKBarrelCollection",
"TPCCollection",
-# "SETCollection",
"OTKBarrelCollection",
"FTDCollection",
"MuonBarrelCollection",
@@ -65,7 +64,7 @@ podioinput = PodioInput("PodioReader", collections=[
## Config ##
vxdhitname = "VXDTrackerHits"
-sithitname = "SITTrackerHits"
+sithitname = "ITKBarrelTrackerHits"
gashitname = "TPCTrackerHits"
sethitname = "OTKBarrelTrackerHits"
setspname = "OTKBarrelSpacePoints"
@@ -86,18 +85,18 @@ digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation"
digiVXD.DigiTool = "SmearDigiTool/VXD"
#digiVXD.OutputLevel = DEBUG
-## SIT ##
-sittool = SmearDigiTool("SIT")
-sittool.ResolutionU = [0.0098]
-sittool.ResolutionV = [0.0433]
-#sittool.OutputLevel = DEBUG
+## ITKBarrel ##
+itkbtool = SmearDigiTool("ITKBarrel")
+itkbtool.ResolutionU = [0.008]
+itkbtool.ResolutionV = [0.040]
+#itkbtool.OutputLevel = DEBUG
-digiSIT = SiTrackerDigiAlg("SITDigi")
-digiSIT.SimTrackHitCollection = "SITCollection"
-digiSIT.TrackerHitCollection = sithitname
-digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
-digiSIT.DigiTool = "SmearDigiTool/SIT"
-#digiSIT.OutputLevel = DEBUG
+digiITKB = SiTrackerDigiAlg("ITKBarrelDigi")
+digiITKB.SimTrackHitCollection = "ITKBarrelCollection"
+digiITKB.TrackerHitCollection = sithitname
+digiITKB.TrackerHitAssociationCollection = "ITKBarrelTrackerHitAssociation"
+digiITKB.DigiTool = "SmearDigiTool/ITKBarrel"
+#digiITKB.OutputLevel = DEBUG
## OTKBarrel ##
otkbtool = SmearDigiTool("OTKBarrel")
@@ -171,8 +170,22 @@ kt110 = KalTestTool("KalTest110")
kt110.Smooth = False
#kt110.OutputLevel = DEBUG
+# Close energy loss
+kt101 = KalTestTool("KalTest101")
+kt101.EnergyLossOn = False
+#kt101.OutputLevel = DEBUG
+
from Configurables import SiliconTrackingAlg
tracking = SiliconTrackingAlg("SiliconTracking")
+tracking.LayerCombinations = [8,7,6, 8,7,5, 8,7,4, 8,7,3, 8,7,2, 8,7,1, 8,7,0,
+ 8,6,5, 8,6,4, 8,6,3, 8,6,2, 8,6,1, 8,6,0,
+ 7,6,5, 7,6,4, 7,6,3, 7,6,2, 7,6,1, 7,6,0,
+ 7,5,3, 7,5,2, 7,5,1, 7,5,0, 7,4,3, 7,4,2, 7,4,1, 7,4,0,
+ 6,5,3, 6,5,2, 6,5,1, 6,5,0, 6,4,3, 6,4,2, 6,4,1, 6,4,0,
+ 6,3,2, 6,3,1, 6,3,0, 6,2,1, 6,2,0, 6,1,0,
+ 5,3,2, 5,3,1, 5,3,0, 5,2,1, 5,2,0, 5,1,0,
+ 4,3,2, 4,3,1, 4,3,0, 4,2,1, 4,2,0, 4,1,0,
+ 3,2,1, 3,2,0, 3,1,0, 2,1,0]
tracking.LayerCombinationsFTD = []
tracking.HeaderCol = "EventHeader"
tracking.VTXHitCollection = vxdhitname
@@ -237,9 +250,8 @@ full.SITRawHits = "NotNeedForPixelSIT"
full.SETRawHits = "NotNeedForPixelSET"
full.FTDRawHits = ftdhitname
full.TPCTracks = "ClupatraTracks" # add standalone TPC track
-full.SiTracks = "SubsetTracks"
+full.SiTracks = "SiTracks"
full.OutputTracks = "CompleteTracks" # default name
-#full.VTXHitToTrackDistance = 5.
full.FTDHitToTrackDistance = 5.
full.SITHitToTrackDistance = 3.
full.SETHitToTrackDistance = 5.
@@ -260,15 +272,16 @@ tof = TofRecAlg("TofRecAlg")
from Configurables import TrackParticleRelationAlg
tpr = TrackParticleRelationAlg("Track2Particle")
tpr.MCParticleCollection = "MCParticle"
-tpr.TrackList = ["CompleteTracks", "ClupatraTracks"]
-tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
+tpr.TrackList = ["CompleteTracks"]
+tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "ITKBarrelTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation"]
#tpr.OutputLevel = DEBUG
+
from Configurables import TrueMuonTagAlg
tmt = TrueMuonTagAlg("TrueMuonTag")
tmt.MCParticleCollection = "MCParticle"
tmt.TrackList = ["CompleteTracks"]
-tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
+tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "ITKBarrelTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%)
tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta)
#tmt.OutputLevel = DEBUG
@@ -282,7 +295,7 @@ out.outputCommands = ["keep *"]
# ApplicationMgr
from Configurables import ApplicationMgr
mgr = ApplicationMgr(
- TopAlg = [podioinput, digiVXD, digiSIT, digiOTKB, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out],
+ TopAlg = [podioinput, digiVXD, digiITKB, digiOTKB, digiFTD, digiTPC, digiMuon, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tof, tmt, out],
EvtSel = 'NONE',
EvtMax = 50,
ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],
diff --git a/Detector/DetCRD/src/Tracker/SiTrackerComposite_v02_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerComposite_v02_geo.cpp
index 35848cb35f835dcdcf91f3bc026fe785af5af771..f2a5a2a791b935d7c274e927cc41f1f2233f6fc9 100644
--- a/Detector/DetCRD/src/Tracker/SiTrackerComposite_v02_geo.cpp
+++ b/Detector/DetCRD/src/Tracker/SiTrackerComposite_v02_geo.cpp
@@ -256,7 +256,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
thisLayer.thicknessSensitive = sensor_thickness;
thisLayer.radiusSupport = flex_radius;
thisLayer.thicknessSupport = flex_thickness;
- //thisLayer.width = flex_width;
+ thisLayer.width = flex_width;
thisLayer.phi0 = phi0 + phi;
thisLayer.rgap = radius;
}
@@ -515,8 +515,8 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
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.distanceSupport = sensitive_radius + support_height / 2.0 - support_thickness / 2.0;
+ topLayer.thicknessSupport = support_thickness / 2.0 + flex_thickness;
topLayer.offsetSupport = -ladder_offset;
topLayer.widthSupport = support_width;
topLayer.zHalfSupport = support_length / 2.0;
@@ -531,7 +531,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
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.thicknessSupport = support_thickness / 2.0 + flex_thickness;
bottomLayer.offsetSupport = -ladder_offset;
bottomLayer.widthSupport = support_width;
bottomLayer.zHalfSupport = support_length / 2.0;
diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v03_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v03_geo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9a268d299b1d6190a6380b4a0bd90f5b3d239eb9
--- /dev/null
+++ b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v03_geo.cpp
@@ -0,0 +1,425 @@
+//====================================================================
+// cepcgeo - CEPC silicon tracker detector models in DD4hep
+//--------------------------------------------------------------------
+// FU Chengdong, 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>
+
+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;
+
+/** Construction of the silicon detector
+ *
+ * ladder staggered to a cirle, like v01 and v02
+ * three supper layer: support, sensitive and flex/power/readout
+ * simple to uniform in each supper layer
+ * \
+ * ___________ \
+ * \ x-y
+ * \
+ * \
+ * sensor: s--side, d--dead
+ * ___________________
+ * | _______s_______ |
+ * | | | |
+ * | |_______d_______| |
+ * | | | |
+ * |s| |s| ----------> z
+ * | | active | |
+ * | | | |
+ * | |_______________| |
+ * |_________s_ _______|
+ * @author FU Chengdong, IHEP, Jan 2025
+ */
+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 printLevel = dd4hep::ERROR;
+ if (x_det.hasAttr(_Unicode(printLevel))) {
+ printLevel = dd4hep::printLevel(x_det.attr<string>(_Unicode(printLevel)));
+ }
+ dd4hep::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 SiTrackerStaggeredLadder_v03 ...");
+
+ dd4hep::rec::ZPlanarData* zPlanarData = new dd4hep::rec::ZPlanarData;
+
+ //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);
+
+ int layer_id = x_layer.attr<int>(_Unicode(layer_id));
+ dd4hep::printout(dd4hep::INFO, "Construct", "layer_id: %02d", layer_id);
+
+ dd4hep::PlacedVolume pv;
+
+ bool dmode = false;
+ if (x_layer.hasAttr(_Unicode(distance))) {
+ dmode = true;
+ }
+ else if (x_layer.hasAttr(_Unicode(radius))) {
+ dmode = false;
+ }
+ else {
+ dd4hep::printout(dd4hep::ERROR, "Construct", "not found parameter distance and radius, at least one");
+ throw runtime_error("dd4hep: lack parameter");
+ }
+
+ std::string layer_name = name + dd4hep::_toString(layer_id, "_layer%02d");
+
+ DetElement ladderDE(layer_name + "_ladder", det_id);
+ DetElement supportDE(layer_name + "_supp", det_id);
+ DetElement sensitiveDE(layer_name + "_sens", det_id);
+ DetElement flexDE(layer_name + "_flex", det_id);
+ DetElement moduleDE(layer_name + "_module", det_id);
+ DetElement sensorDE(layer_name + "_sensor", det_id);
+
+ 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.);
+
+ xml_comp_t x_support(x_layer.child(_Unicode(support)));
+
+ double support_thickness = 0;
+ double support_length = x_support.attr<double>(_Unicode(length));
+ double support_width = x_support.attr<double>(_Unicode(width));
+
+ for (xml_coll_t support_i(x_support,_U(slice)); support_i; ++support_i) {
+ xml_comp_t x_support_slice(support_i);
+ double support_slice_thickness = x_support_slice.attr<double>(_Unicode(thickness));
+ support_thickness += support_slice_thickness;
+ }
+ dd4hep::printout(dd4hep::INFO, "Construct", "support: thickness = %f mm, width = %f mm, length = %f mm",
+ support_thickness/mm, support_width/mm, support_length/mm);
+
+ Material supportMat = theDetector.material(x_support.materialStr());
+ Box supportSolid(support_thickness/2.0, support_width/2.0, support_length/2.0);
+ Volume supportLogical(layer_name + "_SupportLog", supportSolid, supportMat);
+ supportLogical.setVisAttributes(theDetector.visAttributes(x_support.visStr()));
+
+ dd4hep::rec::VolPlane supportSurf(supportLogical, dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Helper),
+ support_thickness/2.0, support_thickness/2.0, u, v, n, o);
+ dd4hep::rec::volSurfaceList(supportDE)->push_back(supportSurf);
+
+ double support_xpos = -support_thickness/2.0;
+ for (xml_coll_t support_i(x_support,_U(slice)); support_i; ++support_i) {
+ xml_comp_t x_support_slice(support_i);
+ double support_slice_thickness = x_support_slice.attr<double>(_Unicode(thickness));
+ double support_slice_width = support_width;
+ double support_slice_length = support_length;
+ double y_offset = 0;
+ double z_offset = 0;
+
+ if (x_support_slice.hasAttr(_Unicode(width))) support_slice_width = x_support_slice.width();//attr<double>(_Unicode(width));
+ if (x_support_slice.hasAttr(_Unicode(length))) support_slice_length = x_support_slice.length();//attr<double>(_Unicode(length));
+ if (x_support_slice.hasAttr(_Unicode(y0))) y_offset = x_support_slice.y0();
+ if (x_support_slice.hasAttr(_Unicode(z0))) z_offset = x_support_slice.z0();
+
+ Material sliceMat = theDetector.material(x_support_slice.materialStr());
+ Box sliceSolid(support_slice_thickness/2.0, support_slice_width/2.0, support_slice_length/2.0);
+ Volume sliceLogical(layer_name + "_SupportSliceLog", sliceSolid, sliceMat);
+ sliceLogical.setVisAttributes(theDetector.visAttributes(x_support_slice.visStr()));
+
+ support_xpos += support_slice_thickness/2.0;
+ pv = supportLogical.placeVolume(sliceLogical, Position(support_xpos, y_offset, z_offset));
+ support_xpos += support_slice_thickness/2.0;
+ }
+
+ xml_comp_t x_sensitive(x_layer.child(_Unicode(sensitive)));
+
+ double sensitive_thickness = x_sensitive.thickness();
+ double sensitive_length = x_sensitive.length();
+ double sensitive_width = x_sensitive.width();
+ double sensitive_center = x_sensitive.hasAttr(_Unicode(gap)) ? x_sensitive.gap() : 0;
+ dd4hep::printout(dd4hep::INFO, "Construct", "sensitive: thickness = %f mm, width = %f mm, length = %f mm",
+ sensitive_thickness/mm, sensitive_width/mm, sensitive_length/mm);
+
+ Material sensitiveMat = theDetector.material(x_sensitive.materialStr());
+ Box sensitiveSolid(sensitive_thickness/2.0, sensitive_width/2.0, sensitive_length/2.0);
+ Volume sensitiveLogical(layer_name + "_SensitiveLog", sensitiveSolid, sensitiveMat);
+ sensitiveLogical.setVisAttributes(theDetector.visAttributes(x_sensitive.visStr()));
+
+ dd4hep::rec::VolPlane sensitiveSurf(sensitiveLogical, dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Helper),
+ sensitive_thickness/2.0, sensitive_thickness/2.0, u, v, n, o);
+ dd4hep::rec::volSurfaceList(sensitiveDE)->push_back(sensitiveSurf);
+
+ xml_comp_t x_module(x_sensitive.child(_Unicode(module)));
+ double module_thickness = sensitive_thickness;
+ double module_length = x_module.length();
+ double module_width = x_module.width();
+ double module_y0 = x_module.hasAttr(_Unicode(y0)) ? x_module.y0() : 0;
+ double module_z0 = x_module.hasAttr(_Unicode(z0)) ? x_module.z0() : 0;
+ dd4hep::printout(dd4hep::INFO, "Construct", "module: thickness = %f mm, width = %f mm, length = %f mm, y0 = %f mm, z0 = %f mm",
+ module_thickness/mm, module_width/mm, module_length/mm, module_y0/mm, module_z0/mm);
+
+ Material moduleMat = theDetector.material(x_module.materialStr());
+ Box moduleSolid(module_thickness/2.0, module_width/2.0, module_length/2.0);
+ Volume moduleLogical(layer_name + "_ModuleLog", moduleSolid, moduleMat);
+ moduleLogical.setVisAttributes(theDetector.visAttributes(x_module.visStr()));
+
+ xml_comp_t x_row(x_module.child(_Unicode(row)));
+ int row_n = x_row.repeat();
+ double row_gap = x_row.gap();
+ double row_side = x_row.attr<double>(_Unicode(up_side));
+ //double row_down_side = x_row.attr<double>(_Unicode(down_side));
+
+ xml_comp_t x_column(x_module.child(_Unicode(column)));
+ int column_n = x_column.repeat();
+ double column_gap = x_column.gap();
+ //double column_side = x_column.attr<double>(_Unicode(left_side));
+ //double column_right_side = x_column.attr<double>(_Unicode(right_side));
+
+ xml_comp_t x_sensor(x_module.child(_Unicode(sensor)));
+ double sensor_thickness = sensitive_thickness;
+ double sensor_length = x_sensor.length();
+ double sensor_width = x_sensor.width();
+ int sensor_n = x_sensor.repeat();
+ double sensor_gap = x_sensor.gap();
+ double sensor_side = x_sensor.attr<double>(_Unicode(side));
+ double sensor_dead = x_sensor.attr<double>(_Unicode(dead));
+ int sensor_row_n = std::floor((module_width-sensor_dead)/sensor_width);
+ int sensor_column_n = sensor_n/sensor_row_n;
+ double sensor_total_length = (sensor_length + sensor_side*2 + sensor_gap) * sensor_column_n - sensor_gap;
+ dd4hep::printout(dd4hep::INFO, "Construct", "sensor: thickness = %f mm, width = %f mm, length = %f mm, total length = %f mm",
+ sensor_thickness/mm, sensor_width/mm, sensor_length/mm, sensor_total_length/mm);
+ dd4hep::printout(dd4hep::INFO, "Construct", " nrow = %d, ncolumn = %d, gap = %f mm, side = %f mm, dead = %f mm",
+ sensor_row_n, sensor_column_n, sensor_gap/mm, sensor_side/mm, sensor_dead/mm);
+
+ Material sensorMat = theDetector.material(x_sensor.materialStr());
+ Box sensorSolid(sensor_thickness/2.0, sensor_width/2.0, sensor_length/2.0);
+ Volume sensorLogical(layer_name + "_SensorLog", sensorSolid, sensorMat);
+ sensorLogical.setVisAttributes(theDetector.visAttributes(x_sensor.visStr()));
+ sensorLogical.setSensitiveDetector(sens);
+ if (x_det.hasAttr(_U(limits))) sensorLogical.setLimitSet(theDetector, x_det.limitsStr());
+
+ dd4hep::rec::VolPlane surf(sensorLogical, dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive),
+ sensor_thickness/2.0, sensor_thickness/2.0, u, v, n, o);
+
+ for (int sensor_id = 0; sensor_id < sensor_n; sensor_id++) {
+ int irow = std::floor(sensor_id/sensor_column_n);
+ int icolumn = sensor_id%sensor_column_n;
+ double ypos = module_y0 - module_width/2.0 + sensor_dead + sensor_side + irow*sensor_gap + (irow + 0.5)*sensor_width;
+ double zpos = module_z0 - sensor_total_length/2.0 + sensor_side + icolumn*(sensor_gap + sensor_side*2) + (icolumn + 0.5)*sensor_length;
+ pv = moduleLogical.placeVolume(sensorLogical, Position(0, ypos, zpos));
+ pv.addPhysVolID("sensor", sensor_id);
+
+ DetElement currentDE = sensorDE.clone(layer_name + dd4hep::_toString(sensor_id, "_sensor%03d"), det_id);
+ currentDE.setPlacement(pv);
+ moduleDE.add(currentDE);
+ dd4hep::rec::volSurfaceList(currentDE)->push_back(surf);
+ }
+
+ double sensitive_active_length = module_length*column_n + sensitive_center + column_gap*column_n;
+ for (int module_id = 0; module_id < row_n*column_n; module_id++) {
+ int irow = std::floor(module_id/column_n);
+ int icolumn = module_id%column_n;
+ double ypos = -sensitive_width/2.0 + row_side + irow*row_gap + (irow + 0.5)*module_width;
+ double zpos = -sensitive_active_length/2.0 + (icolumn + 0.5)*(module_length + column_gap);
+ if (zpos>0) zpos += sensitive_center;
+ pv = sensitiveLogical.placeVolume(moduleLogical, Position(0, ypos, zpos));
+ pv.addPhysVolID("module", module_id);
+
+ DetElement currentDE = moduleDE.clone(layer_name + dd4hep::_toString(module_id, "_module%02d"), det_id);
+ currentDE.setPlacement(pv);
+ sensitiveDE.add(currentDE);
+ }
+
+ xml_comp_t x_flex(x_layer.child(_Unicode(flex)));
+
+ double flex_thickness = 0;
+ double flex_length = x_flex.attr<double>(_Unicode(length));
+ double flex_width = x_flex.attr<double>(_Unicode(width));
+
+ for (xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i) {
+ xml_comp_t x_flex_slice(flex_i);
+ double flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness));
+ flex_thickness += flex_slice_thickness;
+ }
+ dd4hep::printout(dd4hep::INFO, "Construct", "flex: thickness = %f mm, width = %f mm, length = %f mm",
+ flex_thickness/mm, flex_width/mm, flex_length/mm);
+
+ Material flexMat = theDetector.material(x_flex.materialStr());
+ Box flexSolid(flex_thickness/2.0, flex_width/2.0, flex_length/2.0);
+ Volume flexLogical(layer_name + "_FlexLog", flexSolid, flexMat);
+ flexLogical.setVisAttributes(theDetector.visAttributes(x_flex.visStr()));
+
+ dd4hep::rec::VolPlane flexSurf(flexLogical, dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Helper),
+ flex_thickness/2.0, flex_thickness/2.0, u, v, n, o);
+ dd4hep::rec::volSurfaceList(flexDE)->push_back(flexSurf);
+
+ double flex_xpos = -flex_thickness/2.0;
+ for (xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i) {
+ xml_comp_t x_flex_slice(flex_i);
+ double flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness));
+ double flex_slice_width = flex_width;
+ double flex_slice_length = flex_length;
+ double y_offset = 0;
+ double z_offset = 0;
+
+ if (x_flex_slice.hasAttr(_Unicode(width))) flex_slice_width = x_flex_slice.width();//attr<double>(_Unicode(width));
+ if (x_flex_slice.hasAttr(_Unicode(length))) flex_slice_length = x_flex_slice.length();//attr<double>(_Unicode(length));
+ if (x_flex_slice.hasAttr(_Unicode(y0))) y_offset = x_flex_slice.y0();
+ if (x_flex_slice.hasAttr(_Unicode(z0))) z_offset = x_flex_slice.z0();
+
+ Material sliceMat = theDetector.material(x_flex_slice.materialStr());
+ Box sliceSolid(flex_slice_thickness/2.0, flex_slice_width/2.0, flex_slice_length/2.0);
+ Volume sliceLogical(layer_name + "_FlexSliceLog", sliceSolid, sliceMat);
+ sliceLogical.setVisAttributes(theDetector.visAttributes(x_flex_slice.visStr()));
+
+ flex_xpos += flex_slice_thickness/2.0;
+ pv = flexLogical.placeVolume(sliceLogical, Position(flex_xpos, y_offset, z_offset));
+ flex_xpos += flex_slice_thickness/2.0;
+ }
+
+ double ladder_thickness = support_thickness + sensitive_thickness + flex_thickness;
+ double ladder_width = std::max(support_width, std::max(sensitive_width, flex_width));
+ double ladder_length = std::max(support_length, std::max(sensitive_length, flex_length));
+ Box ladderSolid(ladder_thickness/2.0, ladder_width/2.0, ladder_length/2.0);
+ Volume ladderLogical(layer_name + "_LadderLog", ladderSolid, air);
+ ladderLogical.setVisAttributes(theDetector.visAttributes(x_det.visStr()));
+
+ pv = ladderLogical.placeVolume(supportLogical, Position(-(sensitive_thickness + flex_thickness)/2.0, 0, 0));
+ supportDE.setPlacement(pv);
+ ladderDE.add(supportDE);
+
+ pv = ladderLogical.placeVolume(sensitiveLogical, Position( (support_thickness - flex_thickness)/2.0, 0, 0));
+ pv.addPhysVolID("layer", layer_id);
+ sensitiveDE.setPlacement(pv);
+ ladderDE.add(sensitiveDE);
+
+ pv = ladderLogical.placeVolume(flexLogical, Position( (support_thickness + sensitive_thickness)/2.0, 0, 0));
+ flexDE.setPlacement(pv);
+ ladderDE.add(flexDE);
+
+ double ladder_radius = 0;
+ double ladder_distance = 0;
+ double ladder_offset = 0;
+ double ladder_rotate = 0;
+ if (dmode) {
+ ladder_distance = x_layer.attr<double>(_Unicode(distance));
+ ladder_offset = x_layer.attr<double>(_Unicode(offset));
+ ladder_radius = sqrt(ladder_offset*ladder_offset + ladder_distance*ladder_distance);
+ ladder_rotate = -atan(ladder_offset/ladder_distance);
+ }
+ else {
+ double radius = x_layer.attr<double>(_Unicode(radius));
+ //ladder_radius = x_layer.attr<double>(_Unicode(radius));
+ ladder_rotate = x_layer.attr<double>(_Unicode(rotate));
+ //double tanr = fabs(tan(ladder_rotate));
+ //ladder_distance = (tanr*ladder_width/2+sqrt(radius*radius*(tanr*tanr+1)-ladder_width*ladder_width/4))/(tanr*tanr+1);
+ //ladder_distance = ladder_radius*cos(ladder_rotate);
+ ladder_distance = radius*cos(ladder_rotate) - (support_thickness - flex_thickness)/2.0;
+ ladder_radius = ladder_distance/cos(ladder_rotate);
+ //ladder_offset = -ladder_radius*sin(ladder_rotate);
+ ladder_offset = -radius*sin(ladder_rotate);
+ }
+ //int n_sensors_per_ladder = x_layer.attr<int>(_Unicode(n_sensors_per_side));
+ const double ladder_phi0 = x_layer.hasAttr(_Unicode(phi0)) ? x_layer.attr<double>(_Unicode(phi0)) : 0;
+ const int ladder_n = x_layer.attr<int>(_Unicode(n_ladders));
+ const double ladder_dphi = dd4hep::twopi / ladder_n;
+ dd4hep::printout(dd4hep::INFO, "Construct", "ladder: distance = %f mm, offset = %f mm", ladder_distance/mm, ladder_offset/mm);
+ dd4hep::printout(dd4hep::INFO, "Construct", "ladder: center radius = %f mm, phi0 = %f mm", ladder_radius/mm, ladder_phi0);
+ dd4hep::printout(dd4hep::INFO, "Construct", "ladder: n = %d, %d sensors per ladder", ladder_n, row_n*column_n*sensor_n);
+
+ for (int ladder_id = 0; ladder_id < ladder_n; ladder_id++) {
+ DetElement currentDE = ladderDE.clone(layer_name + dd4hep::_toString(ladder_id, "_ladder%03d"), det_id);
+
+ double phi = ladder_phi0 + ladder_dphi*ladder_id;
+ Position pos(ladder_radius*cos(-ladder_rotate+phi), ladder_radius*sin(-ladder_rotate+phi), 0.);
+ Transform3D tr (RotationZYX(phi,0.,0.), pos);
+ pv = envelope.placeVolume(ladderLogical,tr);
+ pv.addPhysVolID("stave", ladder_id);
+ currentDE.setPlacement(pv);
+ det.add(currentDE);
+ }
+
+ // package the reconstruction data
+ dd4hep::rec::ZPlanarData::LayerLayout layer;
+
+ layer.ladderNumber = ladder_n;
+ layer.phi0 = ladder_phi0;
+ layer.sensorsPerLadder = row_n*column_n;
+ layer.lengthSensor = module_length + column_gap;
+ layer.distanceSupport = ladder_distance - ladder_thickness/2.0;
+ layer.thicknessSupport = support_thickness;
+ layer.offsetSupport = ladder_offset;
+ layer.widthSupport = support_width;
+ layer.zHalfSupport = support_length/2.0;
+ layer.distanceSensitive = ladder_distance - ladder_thickness/2.0 + support_thickness;
+ layer.thicknessSensitive = sensitive_thickness;
+ layer.offsetSensitive = ladder_offset;
+ //FIXME: treat widthSensitive as flex_thickness
+ //layer.widthSensitive = module_width;
+ layer.widthSensitive = flex_thickness;
+ layer.zHalfSensitive = sensitive_active_length/2.0;
+
+ zPlanarData->layers.push_back(layer);
+ }
+
+ if (dd4hep::printLevel()<=dd4hep::WARNING) std::cout << (*zPlanarData) << endl;
+ det.addExtension< ZPlanarData >(zPlanarData);
+
+ if (x_det.hasAttr(_U(combineHits))) det.setCombineHits(x_det.attr<bool>(_U(combineHits)),sens);
+
+ dd4hep::printout(dd4hep::INFO, "Construct", "SiTrackerStaggeredLadder_v03 done.");
+ dd4hep::setPrintLevel(oldLevel);
+
+ return det;
+}
+DECLARE_DETELEMENT(SiTrackerStaggeredLadder_v03,create_element)
diff --git a/Detector/DetIdentifier/include/DetIdentifier/CEPCDetectorData.h b/Detector/DetIdentifier/include/DetIdentifier/CEPCDetectorData.h
index e95b62b206ae6fcf1b10fffe32586967663a4423..b466c10cfa3ed9f293727546673aa02d60dc7591 100644
--- a/Detector/DetIdentifier/include/DetIdentifier/CEPCDetectorData.h
+++ b/Detector/DetIdentifier/include/DetIdentifier/CEPCDetectorData.h
@@ -64,6 +64,8 @@ namespace dd4hep {
///
double thicknessSupport;
///
+ double width;
+ ///
double phi0;
///
double rgap;
@@ -88,7 +90,7 @@ namespace dd4hep {
std::vector<CylindricalData::LayerLayout> layers = d.layers;
io << " Layers : " << std::endl
- << " phi0 zHalf radiusSensitive thicknessSensitive radiusSupport thicknessSupport radialgap deltaphi" << std::endl;
+ << " phi0 zHalf width radiusSensitive thicknessSensitive radiusSupport thicknessSupport radialgap deltaphi" << std::endl;
for(unsigned i=0,N=layers.size(); i<N; ++i){
@@ -96,6 +98,7 @@ namespace dd4hep {
io << " " << l.phi0
<< " " << l.zHalf
+ << " " << l.width
<< " " << l.radiusSensitive
<< " " << l.thicknessSensitive
<< " " << l.radiusSupport
@@ -142,13 +145,14 @@ namespace dd4hep {
std::vector<CylindricalData::LayerLayout> layersBent = d.layersBent;
io << " Bent Layers : " << std::endl
- << " phi0 zHalf radiusSensitive thicknessSensitive radiusSupport thicknessSupport radialgap deltaphi" << std::endl;
+ << " phi0 zHalf width 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.width
<< " " << l.radiusSensitive
<< " " << l.thicknessSensitive
<< " " << l.radiusSupport
@@ -187,7 +191,6 @@ namespace dd4hep {
return io;
};
-
/** Simple data structure with key parameters for
* reconstruction of long crystal bar ECAL
*
@@ -201,15 +204,15 @@ namespace dd4hep {
int partNumber;
struct LayerInfo {
- LayerInfo() :
- dlayerNumber(-1),
- slayerNumber(-1),
- barNumber(-1){
- }
- int dlayerNumber;
- int slayerNumber;
- int barNumber;
- };
+ LayerInfo() :
+ dlayerNumber(-1),
+ slayerNumber(-1),
+ barNumber(-1){
+ }
+ int dlayerNumber;
+ int slayerNumber;
+ int barNumber;
+ };
std::vector<LayerInfo> LayerInfos;
};
typedef StructExtension<ECALModuleInfoStruct> ECALModuleInfoData;
@@ -220,7 +223,6 @@ namespace dd4hep {
std::vector<ECALModuleInfoStruct> ModuleInfos;
};
typedef StructExtension<ECALSystemInfoStruct> ECALSystemInfoData;
-
}
}
diff --git a/Digitization/DigiSimple/src/SmearDigiTool.cpp b/Digitization/DigiSimple/src/SmearDigiTool.cpp
index 474b13778ee5d774e327008900e900886abcf51d..64e52e6a28c4d5f35f143d1c3790d4c015e01bc9 100644
--- a/Digitization/DigiSimple/src/SmearDigiTool.cpp
+++ b/Digitization/DigiSimple/src/SmearDigiTool.cpp
@@ -180,10 +180,15 @@ StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHi
<< " outer: " << surface->outerMaterial().name() << endmsg;
if (typeSurface.isPlane() || typeSurface.isCylinder()) {
+ // scale to cov at edge
+ double scale = 1.0;
+
dd4hep::rec::Vector2D localPoint = surface->globalToLocal(oldPos);
+
+ // for planar, same calculation by local is ok, but reduce repeat
if (typeSurface.isCylinder()) {
dd4hep::rec::Vector3D mom_ddrec(mom.x*dd4hep::GeV, mom.y*dd4hep::GeV, mom.z*dd4hep::GeV);
- double length = simhit.getPathLength()*dd4hep::mm;
+ double length = simhit.getPathLength()*dd4hep::mm;
dd4hep::rec::Vector3D pre = oldPos - (0.5*length)*mom_ddrec.unit();
dd4hep::rec::Vector3D post = oldPos + (0.5*length)*mom_ddrec.unit();
dd4hep::rec::Vector2D localPre = surface->globalToLocal(pre);
@@ -192,14 +197,36 @@ StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHi
debug() << "pre: (" << pre.x() << " " << pre.y() << " " << pre.z() << " ) local (" << localPre.u() << ", " << localPre.v() << " ) "
<< "post: (" << post.x() << " " << post.y() << " " << post.z() << " ) local (" << localPost.u() << ", " << localPost.v() << " ) " << endmsg;
}
- dd4hep::rec::Vector3D local3D(localPoint.u(), localPoint.v(), 0);
+ //dd4hep::rec::Vector3D local3D(localPoint.u(), localPoint.v(), 0);
// A small check, if the hit is in the boundaries:
if (!surface->insideBounds(oldPos)) {
- warning() << " global: (" << oldPos.x() << " " << oldPos.y() << " " << oldPos.z()
- << ") local: (" << localPoint.u() << ", " << localPoint.v() << " )"
- << " step: " << simhit.getPathLength()
- << " is not within boundaries. Hit is skipped." << endmsg;
- return StatusCode::SUCCESS;
+ double dSToHit = surface->distance(oldPos);
+ debug() << " global: (" << oldPos.x()/dd4hep::mm << " " << oldPos.y()/dd4hep::mm << " " << oldPos.z()/dd4hep::mm
+ << ") local: (" << localPoint.u()/dd4hep::mm << ", " << localPoint.v()/dd4hep::mm << " )"
+ << " distance: " << dSToHit/dd4hep::mm
+ << " is not within boundaries." << endmsg;
+
+ // FIXME: change position to path at center plane? or enlarge cov? other tracker
+ if (system == CEPCConf::DetID::OTKBarrel || system == CEPCConf::DetID::OTKEndcap) {
+#if 0
+ dd4hep::rec::Vector3D global3D = oldPos;
+ dd4hep::rec::Vector3D mom_ddrec(mom.x*dd4hep::GeV, mom.y*dd4hep::GeV, mom.z*dd4hep::GeV);
+ double length = simhit.getPathLength()*dd4hep::mm;
+ dd4hep::rec::Vector3D pre = oldPos - (0.5*length)*mom_ddrec.unit();
+ dd4hep::rec::Vector3D post = oldPos + (0.5*length)*mom_ddrec.unit();
+ double dSToPre = surface->distance(pre);
+ double dPreToHit = dSToPre - dSToHit;
+ global3D = pre + (length/dPreToHit*dSToPre)*mom_ddrec.unit();
+ localPoint = surface->globalToLocal(global3D);
+ scale = 1;
+
+ debug() << " global: (" << global3D.x()/dd4hep::mm << " " << global3D.y()/dd4hep::mm << " " << global3D.z()/dd4hep::mm
+ << ") local: (" << localPoint.u()/dd4hep::mm << ", " << localPoint.v()/dd4hep::mm << ", 0 )"
+ << " distance: " << surface->distance(global3D)/dd4hep::mm
+ << " is not within boundaries." << endmsg;
+#endif
+ }
+ //else return StatusCode::SUCCESS;
}
dd4hep::rec::Vector3D globalPointSmeared;//CLHEP::Hep3Vector globalPoint(pos[0],pos[1],pos[2]);
dd4hep::rec::Vector2D localPointSmeared;
@@ -265,19 +292,19 @@ StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHi
std::array<float, 6> cov;
cov[0] = u_direction[0];
cov[1] = u_direction[1];
- cov[2] = resU;
+ cov[2] = resU*scale;
cov[3] = v_direction[0];
cov[4] = v_direction[1];
- cov[5] = resV;
+ cov[5] = resV*scale;
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]));
+ outhit.setCovMatrix(CEPC::ConvertToCovXYZ(resU*scale, u_direction[0], u_direction[1], resV*scale, 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});
+ outhit.setCovMatrix(std::array<float, 6>{resU*resU*scale*scale/2., 0, resU*resU*scale*scale/2, 0, 0, resV*resV*scale*scale});
type.set(CEPCConf::TrkHitTypeBit::CYLINDER);
}
diff --git a/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp b/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp
index b2972273323733713eae3e3eecd54f6ea72521dd..6df266780ec1b847ad1ccf62600cf27d85d859e4 100644
--- a/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp
+++ b/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.cpp
@@ -4,6 +4,9 @@
#include "TrackSystemSvc/MarlinTrkUtils.h"
#include "TrackSystemSvc/IMarlinTrack.h"
#include "DetInterface/IGeomSvc.h"
+#include "DetIdentifier/CEPCConf.h"
+
+#include "DataHelper/Navigation.h"
#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
@@ -11,6 +14,15 @@
#include "edm4hep/TrackerHit.h"
#include "edm4hep/TrackState.h"
#include "edm4hep/MutableTrack.h"
+//#if __has_include("edm4hep/EDM4hepVersion.h")
+//#include "edm4hep/EDM4hepVersion.h"
+//#else
+// Copy the necessary parts from the header above to make whatever we need to work here
+//#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch)))
+// v00-09 is the last version without the capitalization change of the track vector members
+//#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0)
+//#endif
+#include "podio/podioVersion.h"
DECLARE_COMPONENT(KalTestTool)
@@ -69,7 +81,7 @@ int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector<edm4hep::TrackerHi
debug() << "start..." << endmsg;
std::shared_ptr<MarlinTrk::IMarlinTrack> marlinTrack(m_factoryMarlinTrk->createTrack());
debug() << "created MarlinKalTestTrack" << endmsg;
- int status = MarlinTrk::createFinalisedLCIOTrack(marlinTrack.get(), trackHits, &track, backward, covMatrix, m_magneticField, maxChi2perHit);
+ int status = this->createFinalisedTrack(marlinTrack.get(), trackHits, &track, backward, covMatrix, m_magneticField, maxChi2perHit);
marlinTrack->getHitsInFit(m_hitsInFit);
marlinTrack->getOutliers(m_outliers);
@@ -91,7 +103,7 @@ int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector<edm4hep::TrackerHi
debug() << "start..." << endmsg;
std::shared_ptr<MarlinTrk::IMarlinTrack> marlinTrack(m_factoryMarlinTrk->createTrack());
debug() << "created MarlinKalTestTrack" << endmsg;
- int status = MarlinTrk::createFinalisedLCIOTrack(marlinTrack.get(), trackHits, &track, backward, &trackState, m_magneticField, maxChi2perHit);
+ int status = this->createFinalisedTrack(marlinTrack.get(), trackHits, &track, backward, &trackState, m_magneticField, maxChi2perHit);
marlinTrack->getHitsInFit(m_hitsInFit);
marlinTrack->getOutliers(m_outliers);
@@ -101,3 +113,409 @@ int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector<edm4hep::TrackerHi
error() << "Don't support the Fitter " << m_fitterName << endmsg;
return 0;
}
+
+// from MarlinTrk
+int KalTestTool::createFinalisedTrack(MarlinTrk::IMarlinTrack* marlinTrk, std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::MutableTrack* track, bool fit_backwards,
+ const decltype(edm4hep::TrackState::covMatrix)& initial_cov_for_prefit, float bfield_z, double maxChi2Increment) {
+ ///////////////////////////////////////////////////////
+ // check inputs
+ ///////////////////////////////////////////////////////
+ if ( hit_list.empty() ) return MarlinTrk::IMarlinTrack::bad_intputs;
+
+ if( track == 0 ){
+ throw std::runtime_error( "MarlinTrk::finaliseLCIOTrack: TrackImpl == NULL " ) ;
+ }
+
+ int return_error = 0;
+ ///////////////////////////////////////////////////////
+ // produce prefit parameters
+ ///////////////////////////////////////////////////////
+
+ edm4hep::TrackState pre_fit ;
+
+ //std::cout << "debug:=====================before createPrefit" << std::endl;
+ return_error = MarlinTrk::createPrefit(hit_list, &pre_fit, bfield_z, fit_backwards);
+ //std::cout << "debug:=====================after createPrefit return code=" << return_error << std::endl;
+ pre_fit.covMatrix = initial_cov_for_prefit;
+
+ ///////////////////////////////////////////////////////
+ // use prefit parameters to produce Finalised track
+ ///////////////////////////////////////////////////////
+
+ if( return_error == 0 ) {
+
+ return_error = createFinalisedTrack( marlinTrk, hit_list, track, fit_backwards, &pre_fit, bfield_z, maxChi2Increment);
+
+ }
+ else {
+ warning() << "createFinalisedLCIOTrack : Prefit failed error = " << return_error << endmsg;
+ }
+ return return_error;
+}
+
+// from MarlinTrk
+int KalTestTool::createFinalisedTrack(MarlinTrk::IMarlinTrack* marlinTrk, std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::MutableTrack* track, bool fit_backwards,
+ edm4hep::TrackState* pre_fit, float bfield_z, double maxChi2Increment) {
+
+ ///////////////////////////////////////////////////////
+ // check inputs
+ ///////////////////////////////////////////////////////
+ if ( hit_list.empty() ) return MarlinTrk::IMarlinTrack::bad_intputs;
+
+ if( track == 0 ){
+ throw std::runtime_error("MarlinTrk::finaliseLCIOTrack: TrackImpl == NULL ");
+ }
+
+ if( pre_fit == 0 ){
+ throw std::runtime_error("MarlinTrk::finaliseLCIOTrack: TrackStateImpl == NULL ");
+ }
+
+
+ int fit_status = MarlinTrk::createFit(hit_list, marlinTrk, pre_fit, bfield_z, fit_backwards, maxChi2Increment);
+
+ if (fit_status != MarlinTrk::IMarlinTrack::success) {
+ debug() << "createFinalisedTrack fit failed: fit_status = " << fit_status << endmsg;
+ return fit_status;
+ }
+
+ int error = finaliseTrack(marlinTrk, track, hit_list, fit_backwards);
+ debug() << "finaliseTrack. status = " << error << endmsg;
+
+ return error;
+}
+
+// from MarlinTrk
+int KalTestTool::finaliseTrack(MarlinTrk::IMarlinTrack* marlintrk, edm4hep::MutableTrack* track, std::vector<edm4hep::TrackerHit>& hit_list, bool fit_backwards,
+ edm4hep::TrackState* atLastHit, edm4hep::TrackState* atCaloFace) {
+
+ ///////////////////////////////////////////////////////
+ // check inputs
+ ///////////////////////////////////////////////////////
+ if (marlintrk == 0) {
+ throw std::runtime_error("MarlinTrk::finaliseLCIOTrack: IMarlinTrack == NULL ");
+ }
+
+ if (track == 0) {
+ throw std::runtime_error("MarlinTrk::finaliseLCIOTrack: TrackImpl == NULL ");
+ }
+
+ if (atCaloFace && atLastHit == 0) {
+ throw std::runtime_error("MarlinTrk::finaliseLCIOTrack: atLastHit == NULL ");
+ }
+
+ if (atLastHit && atCaloFace == 0) {
+ throw std::runtime_error("MarlinTrk::finaliseLCIOTrack: atCaloFace == NULL ");
+ }
+
+ ///////////////////////////////////////////////////////
+ // error to return if any
+ ///////////////////////////////////////////////////////
+ int return_error = 0;
+
+ int ndf = 0;
+ double chi2 = -DBL_MAX;
+
+ /////////////////////////////////////////////////////////////
+ // First check NDF to see if it make any sense to continue.
+ // The track will be dropped if the NDF is less than 0
+ /////////////////////////////////////////////////////////////
+
+ return_error = marlintrk->getNDF(ndf);
+
+ if (return_error != MarlinTrk::IMarlinTrack::success) {
+ debug() << "getNDF returns " << return_error << endmsg;
+ return return_error;
+ }
+ else if (ndf < 0) {
+ debug() << "number of degrees of freedom less than 0 track dropped : NDF = " << ndf << endmsg;
+ return MarlinTrk::IMarlinTrack::error;
+ }
+ else {
+ debug() << "NDF = " << ndf << endmsg;
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // get the list of hits used in the fit
+ // add these to the track, add spacepoints as long as at least on strip hit is used.
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ std::vector<std::pair<edm4hep::TrackerHit, double> > hits_in_fit;
+ std::vector<std::pair<edm4hep::TrackerHit, double> > outliers;
+ std::vector<edm4hep::TrackerHit> used_hits;
+
+ hits_in_fit.reserve(300);
+ outliers.reserve(300);
+
+ marlintrk->getHitsInFit(hits_in_fit);
+ marlintrk->getOutliers(outliers);
+
+ ///////////////////////////////////////////////
+ // now loop over the hits provided for fitting
+ // we do this so that the hits are added in the
+ // order in which they have been fitted
+ ///////////////////////////////////////////////
+
+ for ( unsigned ihit = 0; ihit < hit_list.size(); ++ihit) {
+
+ edm4hep::TrackerHit trkHit = hit_list[ihit];
+
+ std::bitset<32> type = trkHit.getType();
+ if (type.test(CEPCConf::TrkHitTypeBit::COMPOSITE_SPACEPOINT)) { //it is a composite spacepoint
+ //std::cout << "Error: space point is not still valid! pelease wait updating..." <<std::endl;
+ //exit(1);
+ // get strip hits
+ int nRawHit = trkHit.rawHits_size();
+ for( unsigned k=0; k< nRawHit; k++ ){
+ edm4hep::TrackerHit rawHit = Navigation::Instance()->GetTrackerHit(trkHit.getRawHits(k));
+ bool is_outlier = false;
+ // here we loop over outliers as this will be faster than looping over the used hits
+ for ( unsigned ohit = 0; ohit < outliers.size(); ++ohit) {
+ if ( rawHit.id() == outliers[ohit].first.id() ) {
+ is_outlier = true;
+ break; // break out of loop over outliers
+ }
+ }
+ if (is_outlier == false) {
+ used_hits.push_back(hit_list[ihit]);
+ track->addToTrackerHits(used_hits.back());
+ break; // break out of loop over rawObjects
+ }
+ }
+ }
+ else {
+ bool is_outlier = false;
+ // here we loop over outliers as this will be faster than looping over the used hits
+ for ( unsigned ohit = 0; ohit < outliers.size(); ++ohit) {
+ if ( trkHit == outliers[ohit].first ) {
+ is_outlier = true;
+ break; // break out of loop over outliers
+ }
+ }
+ if (is_outlier == false) {
+ used_hits.push_back(hit_list[ihit]);
+ track->addToTrackerHits(used_hits.back());
+ }
+ }
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ // We now need to find out at which point the fit is constrained
+ // and therefore be able to provide well formed (pos. def.) cov. matrices
+ ///////////////////////////////////////////////////////////////////////////
+
+ ///////////////////////////////////////////////////////
+ // first hit
+ ///////////////////////////////////////////////////////
+
+ edm4hep::TrackState* trkStateAtFirstHit = new edm4hep::TrackState() ;
+ edm4hep::TrackerHit firstHit = (fit_backwards == MarlinTrk::IMarlinTrack::backward) ? hits_in_fit.back().first : hits_in_fit.front().first;
+
+ ///////////////////////////////////////////////////////
+ // last hit
+ ///////////////////////////////////////////////////////
+
+ edm4hep::TrackState* trkStateAtLastHit = new edm4hep::TrackState() ;
+
+ edm4hep::TrackerHit lastHit = (fit_backwards == MarlinTrk::IMarlinTrack::backward) ? hits_in_fit.front().first : hits_in_fit.back().first;
+
+#if PODIO_BUILD_VERSION < PODIO_VERSION(0, 17, 4)
+ edm4hep::TrackerHit last_constrained_hit(0);
+#else
+ auto last_constrained_hit = edm4hep::TrackerHit::makeEmpty();
+#endif
+
+ marlintrk->getTrackerHitAtPositiveNDF(last_constrained_hit);
+
+ return_error = marlintrk->smooth(lastHit);
+
+ if (return_error != MarlinTrk::IMarlinTrack::success) {
+ debug() << "return_code for smoothing to " << lastHit << " = " << return_error << " NDF = " << ndf << endmsg;
+ delete trkStateAtFirstHit;
+ delete trkStateAtLastHit;
+ return return_error ;
+ }
+
+ ///////////////////////////////////////////////////////
+ // first create trackstate at IP
+ ///////////////////////////////////////////////////////
+ const edm4hep::Vector3d point; // nominal IP
+
+ edm4hep::TrackState* trkStateIP = new edm4hep::TrackState();
+
+ ///////////////////////////////////////////////////////
+ // make sure that the track state can be propagated to the IP
+ ///////////////////////////////////////////////////////
+
+ //FIXME: AidaTT not used now, if to use AidaTT, these code need to update and apply
+ //MarlinTrk::IMarlinTrkSystem* trksystem = MarlinTrk::Factory::getCurrentMarlinTrkSystem() ;
+ //bool usingAidaTT = ( trksystem->name() == "AidaTT" ) ;
+ bool usingAidaTT = false;
+ if (fit_backwards == MarlinTrk::IMarlinTrack::backward || usingAidaTT ) {
+ return_error = marlintrk->propagate(point, firstHit, *trkStateIP, chi2, ndf ) ;
+ }
+ else {
+ // if we fitted forward, we start from the last_constrained hit
+ // and then add the last inner hits with a Kalman step ...
+
+ // create a temporary IMarlinTrack
+ auto mTrk = std::shared_ptr<MarlinTrk::IMarlinTrack>( m_factoryMarlinTrk->createTrack() );
+
+ edm4hep::TrackState ts;
+
+ auto hI = hits_in_fit.rbegin();
+
+ double chi2Tmp = 0;
+ int ndfTmp = 0;
+ //return_error = marlintrk->getTrackState( last_constrained_hit, ts, chi2, ndf);
+ return_error = marlintrk->getTrackState( lastHit, ts, chi2, ndf);
+
+ debug() << "-- TrackState at last constrained hit : " << ts << endmsg;
+
+ //need to add a dummy hit to the track
+ //mTrk->addHit(last_constrained_hit);
+ mTrk->addHit(lastHit);
+
+ double _bfield = m_magneticField;
+ // fixme: the implementation for DDKalTest does no longer need this value but the IMarlinTrk interface is not yet changed
+ mTrk->initialise(ts, _bfield, fit_backwards);
+
+ //while (hI->first.id() != last_constrained_hit.id()) {
+ while (hI->first.id() != lastHit.id()) {
+ debug() << "-- hit in reverse_iterator : " << hI->first.getCellID() << " " << hI->first.getPosition() << endmsg;
+ ++hI;
+ }
+
+ ++hI;
+
+ while (hI != hits_in_fit.rend()) {
+ auto hit = (*hI).first;
+
+ double deltaChi;
+ double maxChi2Increment = DBL_MAX; // not apply chisquare cut, since hits_in_fit have past in filter
+
+ int addHit = mTrk->addAndFit(hit, deltaChi, maxChi2Increment);
+
+ debug() << "-- hit " << hit << " added : " << MarlinTrk::errorCode(addHit)
+ << " deltaChi2: " << deltaChi << endmsg;
+
+ if (addHit != MarlinTrk::IMarlinTrack::success) {
+ debug() << "-- could not add inner hit to track !!! " << maxChi2Increment << endmsg;
+ }
+
+ ++hI;
+ }//------------------------------------
+
+ debug() << "-- temporary kaltest track for track state at the IP: " << mTrk->toString() << endmsg;
+
+ // now propagate the temporary track to the IP
+ return_error = mTrk->propagate(point, firstHit, *trkStateIP, chi2Tmp, ndfTmp);
+
+ debug() << "-- propagated temporary track fromfirst hit to IP : " << (*trkStateIP) << endmsg;
+ //FIXME: if forward, better by add the last inner hits with a Kalman step from the last_constrained hit
+ //return_error = marlintrk->propagate(point, firstHit, *trkStateIP, chi2, ndf ) ;
+ }
+ if (return_error != MarlinTrk::IMarlinTrack::success) {
+ debug() << "-- return_code for propagation = " << MarlinTrk::errorCode(return_error) << " NDF = " << ndf << std::endl;
+ delete trkStateIP;
+ delete trkStateAtFirstHit;
+ delete trkStateAtLastHit;
+
+ return return_error ;
+ }
+
+ trkStateIP->location = edm4hep::TrackState::AtIP;
+ track->addToTrackStates(*trkStateIP);
+ track->setChi2(chi2);
+ track->setNdf(ndf);
+
+ ///////////////////////////////////////////////////////
+ // set the track states at the first and last hits
+ ///////////////////////////////////////////////////////
+
+ ///////////////////////////////////////////////////////
+ // @ first hit
+ ///////////////////////////////////////////////////////
+
+ return_error = marlintrk->getTrackState(firstHit, *trkStateAtFirstHit, chi2, ndf ) ;
+
+ if ( return_error == 0 ) {
+ trkStateAtFirstHit->location = edm4hep::TrackState::AtFirstHit;
+ track->addToTrackStates(*trkStateAtFirstHit);
+ debug() << ">>>> create TrackState AtFirstHit" << (*trkStateAtFirstHit) << endmsg;
+ }
+ else {
+ warning() << ">>>> could not get TrackState at First Hit " << firstHit.getCellID() << " " << firstHit.getPosition() << endmsg;
+ }
+
+ double r_first = firstHit.getPosition()[0]*firstHit.getPosition()[0] + firstHit.getPosition()[1]*firstHit.getPosition()[1];
+
+ track->setRadiusOfInnermostHit(sqrt(r_first));
+
+ if ( atLastHit == 0 && atCaloFace == 0 ) {
+
+ ///////////////////////////////////////////////////////
+ // @ last hit
+ ///////////////////////////////////////////////////////
+
+ debug() << ">>>> create TrackState AtLastHit : using from trkhit "
+ << last_constrained_hit.getCellID() << " " << last_constrained_hit.getPosition() << endmsg;
+
+ edm4hep::Vector3d last_hit_pos(lastHit.getPosition());
+
+ return_error = marlintrk->propagate(last_hit_pos, last_constrained_hit, *trkStateAtLastHit, chi2, ndf);
+
+ if ( return_error == 0 ) {
+ trkStateAtLastHit->location = edm4hep::TrackState::AtLastHit;
+ track->addToTrackStates(*trkStateAtLastHit);
+
+ debug() << ">>>> createTrackStateAtLastHit OK: " << (*trkStateAtLastHit) << endmsg;
+ }
+ else {
+ error() << ">>>> could not get TrackState at Last Hit " << lastHit.getCellID() << " " << lastHit.getPosition()
+ << " from last_constrained_hit " << last_constrained_hit.getCellID() << " " << last_constrained_hit.getPosition() << endmsg;
+ //delete trkStateAtLastHit;
+ }
+
+ ///////////////////////////////////////////////////////
+ // set the track state at Calo Face
+ ///////////////////////////////////////////////////////
+
+ edm4hep::TrackState trkStateCalo;
+ bool tanL_is_positive = trkStateIP->tanLambda >0 ;
+
+ return_error = MarlinTrk::createTrackStateAtCaloFace(marlintrk, &trkStateCalo, last_constrained_hit, tanL_is_positive);
+
+ if (return_error == 0) {
+ trkStateCalo.location = edm4hep::TrackState::AtCalorimeter;
+ track->addToTrackStates(trkStateCalo);
+
+ debug() << ">>>> createTrackStateAtCaloFace OK: " << trkStateCalo << endmsg;
+ }
+ else {
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << ">>>> could not get TrackState at Calo Face " << endmsg;
+ if (last_constrained_hit.isAvailable()) {
+ auto pos = last_constrained_hit.getPosition();
+ debug() << ">>>> last_constrained_hit = " << pos.x << "," << pos.y << "," << pos.z << endmsg;
+ }
+ else {
+ debug() << ">>>> last_constrained_hit not Available" << endmsg;
+ }
+ }
+ //delete trkStateCalo;
+ }
+ }
+ else {
+ track->addToTrackStates(*atLastHit);
+ track->addToTrackStates(*atCaloFace);
+ //delete trkStateAtLastHit;
+ }
+
+ if(trkStateAtFirstHit) delete trkStateAtFirstHit;
+ if(trkStateAtLastHit) delete trkStateAtLastHit;
+ if(trkStateIP) delete trkStateIP;
+ ///////////////////////////////////////////////////////
+ // done
+ ///////////////////////////////////////////////////////
+ return return_error;
+}
diff --git a/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.h b/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.h
index 0cce58bc46c56d8bc15c0a4670205e58c27aaf38..5aff0859fe10d9b5b66befffdc94cef210095f9c 100644
--- a/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.h
+++ b/Reconstruction/RecTrkGlobal/src/FitterTool/KalTestTool.h
@@ -6,7 +6,8 @@
#include <vector>
namespace MarlinTrk {
- class IMarlinTrkSystem ;
+ class IMarlinTrkSystem;
+ class IMarlinTrack;
}
class KalTestTool : public extends<AlgTool, ITrackFitterTool> {
@@ -26,6 +27,12 @@ class KalTestTool : public extends<AlgTool, ITrackFitterTool> {
std::vector<std::pair<edm4hep::TrackerHit, double> >& GetOutliers() override {return m_outliers;};
void Clear() override {m_hitsInFit.clear(); m_outliers.clear();};
+ int createFinalisedTrack(MarlinTrk::IMarlinTrack* marlinTrk, std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::MutableTrack* track, bool fit_backwards,
+ const decltype(edm4hep::TrackState::covMatrix)& initial_cov_for_prefit, float bfield_z, double maxChi2Increment);
+ int createFinalisedTrack(MarlinTrk::IMarlinTrack* marlinTrk, std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::MutableTrack* track, bool fit_backwards,
+ edm4hep::TrackState* pre_fit, float bfield_z, double maxChi2Increment);
+ int finaliseTrack(MarlinTrk::IMarlinTrack* marlintrk, edm4hep::MutableTrack* track, std::vector<edm4hep::TrackerHit>& hit_list, bool fit_backwards,
+ edm4hep::TrackState* atLastHit=0, edm4hep::TrackState* atCaloFace=0);
private:
Gaudi::Property<std::string> m_fitterName{this, "Fitter", "KalTest"};
Gaudi::Property<bool> m_useQMS{this, "MSOn", true};
diff --git a/Reconstruction/RecTrkGlobal/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/RecTrkGlobal/src/FullLDCTracking/FullLDCTrackingAlg.cpp
index 61160c99816dcd450c9ed0af757dfb09f979a912..001971edfbba00611f1319dd1237153fd0ac48b1 100755
--- a/Reconstruction/RecTrkGlobal/src/FullLDCTracking/FullLDCTrackingAlg.cpp
+++ b/Reconstruction/RecTrkGlobal/src/FullLDCTracking/FullLDCTrackingAlg.cpp
@@ -68,6 +68,9 @@
#include <TStopwatch.h>
#include "TMath.h"
+#include "streamlog/streamlog.h"
+//streamlog::logstream out;
+
typedef std::vector<edm4hep::TrackerHit> TrackerHitVec;
using namespace edm4hep ;
@@ -149,6 +152,9 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc
*/
// Output track collection
declareProperty("OutputTracks", _OutputTrackColHdl, "Full LDC track collection name");
+
+ //out.init(std::cout, name);
+ //out.addLevelName<streamlog::DEBUG>();
}
@@ -217,8 +223,6 @@ void FullLDCTrackingAlg::processRunHeader( LCRunHeader* run) {
*/
StatusCode FullLDCTrackingAlg::execute() {
-
- // debug() << endmsg;
info() << "FullLDCTrackingAlg -> run = " << 0/*evt->getRunNumber()*/ << " event = " << _nEvt << endmsg;
auto stopwatch = TStopwatch();
@@ -420,6 +424,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
int fit_count = 0;
bool sort_by_r = false;
bool use_ts_initial = false;
+ double maxChi2PerHit = _maxChi2PerHit;
fitstart:
@@ -473,10 +478,10 @@ fitstart:
debug() << "fit direction: " << ((fit_backwards==IMarlinTrack::backward) ? "backward" : "forward") << endmsg;
//debug() << "MaxChi2PerHit: " << _maxChi2PerHit << endmsg;
if (use_ts_initial) {
- error_code = m_fitTool->Fit(track, trkHits, ts_initial, _maxChi2PerHit, fit_backwards);
+ error_code = m_fitTool->Fit(track, trkHits, ts_initial, maxChi2PerHit, fit_backwards);
}
else {
- error_code = m_fitTool->Fit(track, trkHits, covMatrix, _maxChi2PerHit, fit_backwards);
+ error_code = m_fitTool->Fit(track, trkHits, covMatrix, maxChi2PerHit, fit_backwards);
}
} catch (...) {
// delete track;
@@ -494,6 +499,14 @@ fitstart:
}
#endif
+ std::vector<std::pair<edm4hep::TrackerHit , double> > hits_in_fit ;
+ std::vector<std::pair<edm4hep::TrackerHit , double> > outliers ;
+
+ UTIL::BitField64 cellID_encoder( lcio::ILDCellID0::encoder_string );
+
+ hits_in_fit = m_fitTool->GetHitsInFit();
+ outliers = m_fitTool->GetOutliers();
+
if (error_code != IMarlinTrack::success) {
debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit failed with error code " << error_code << " track dropped. Number of hits = "<< trkHits.size() << endmsg;
m_fitTool->Clear();
@@ -504,25 +517,23 @@ fitstart:
// TODO: to pick up old tracks
continue ;
}
-
- std::vector<std::pair<edm4hep::TrackerHit , double> > hits_in_fit ;
- std::vector<std::pair<edm4hep::TrackerHit , double> > outliers ;
+ else {
+ // FIXME: better to loose cut?
+ if (fit_count < 3 && float(outliers.size())/float(hits_in_fit.size()) > 0.2) {
+ maxChi2PerHit *= 2;
+ m_fitTool->Clear();
+ goto fitstart;
+ }
+ }
+
std::vector<edm4hep::TrackerHit> all_hits;
all_hits.reserve(hits_in_fit.size());
-
- //marlinTrk->getHitsInFit(hits_in_fit);
- hits_in_fit = m_fitTool->GetHitsInFit();
for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) {
all_hits.push_back(hits_in_fit[ihit].first);
}
-
- UTIL::BitField64 cellID_encoder( lcio::ILDCellID0::encoder_string ) ;
-
+
MarlinTrk::addHitNumbersToTrack(&track, all_hits, true, cellID_encoder);
-
- //marlinTrk->getOutliers(outliers);
- outliers = m_fitTool->GetOutliers();
for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) {
all_hits.push_back(outliers[ihit].first);
@@ -1373,7 +1384,7 @@ void FullLDCTrackingAlg::prepareVectors() {
trackExt->setD0(getD0(tpcTrack));
trackExt->setZ0(getZ0(tpcTrack));
- float cov[15];
+ //float cov[15];
//float param[5];
// float reso[4];
//param[0] = getOmega(tpcTrack);
@@ -1384,11 +1395,13 @@ void FullLDCTrackingAlg::prepareVectors() {
auto Cov = getCovMatrix(tpcTrack);
int NC = int(Cov.size());
+
+ std::unique_ptr<float[]> cov(new float[NC]);
for (int ic=0;ic<NC;ic++) {
cov[ic] = Cov[ic];
}
- trackExt->setCovMatrix(cov);
+ trackExt->setCovMatrix(cov.get());
trackExt->setNDF(tpcTrack.getNdf());
trackExt->setChi2(tpcTrack.getChi2());
for (int iHit=0;iHit<nHits;++iHit) {
@@ -1442,7 +1455,7 @@ void FullLDCTrackingAlg::prepareVectors() {
trackExt->setPhi(getPhi(siTrack));
trackExt->setD0(getD0(siTrack));
trackExt->setZ0(getZ0(siTrack));
- float cov[15];
+ //float cov[15];
//float param[5];
//param[0] = getOmega(siTrack);
@@ -1453,10 +1466,12 @@ void FullLDCTrackingAlg::prepareVectors() {
auto Cov = getCovMatrix(siTrack);
int NC = int(Cov.size());
+
+ std::unique_ptr<float[]> cov(new float[NC]);
for (int ic=0;ic<NC;ic++) {
cov[ic] = Cov[ic];
}
- trackExt->setCovMatrix(cov);
+ trackExt->setCovMatrix(cov.get());
trackExt->setNDF(siTrack.getNdf());
trackExt->setChi2(siTrack.getChi2());
char strg[200];
@@ -1646,7 +1661,6 @@ void FullLDCTrackingAlg::MergeTPCandSiTracks() {
if (combinedTrack != NULL) {
-
_allCombinedTracks.push_back( combinedTrack );
_candidateCombinedTracks.insert(tpcTrackExt);
_candidateCombinedTracks.insert(siTrackExt);
@@ -1843,12 +1857,13 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
int errorCode = IMarlinTrack::success;
try{
- pre_fit = getTrackStateAt(tpcTrack->getTrack(), edm4hep::TrackState::AtLastHit);
+ pre_fit = getTrackStateAt(siTrack->getTrack(), edm4hep::TrackState::AtFirstHit);
+ debug() << pre_fit << endmsg;
}
catch(std::runtime_error& e){
error() << e.what() << " should be checked (TPC track) " << tpcTrack << endmsg;
try{
- pre_fit = getTrackStateAt(siTrack->getTrack(), edm4hep::TrackState::AtLastHit);
+ pre_fit = getTrackStateAt(tpcTrack->getTrack(), edm4hep::TrackState::AtLastHit);
}
catch(std::runtime_error& e){
error() << e.what() << " should be checked (Si track)" << endmsg;
@@ -1870,17 +1885,27 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac
pre_fit.covMatrix = covMatrix;
- errorCode = MarlinTrk::createFit( trkHits, &marlin_trk, &pre_fit, _bField, IMarlinTrack::backward , _maxChi2PerHit );
+ errorCode = MarlinTrk::createFit( trkHits, &marlin_trk, &pre_fit, _bField, !IMarlinTrack::backward , _maxChi2PerHit );
if ( errorCode != IMarlinTrack::success ) {
debug() << "FullLDCTrackingAlg::CombineTracks: creation of fit fails with error " << errorCode << endmsg;
return 0;
}
debug() << "createFit finished" << endmsg;
+
+ std::vector<std::pair<edm4hep::TrackerHit, double> > hitsInFit;
+ marlin_trk.getHitsInFit(hitsInFit);
+ auto firstHit = hitsInFit.begin()->first;
+
+ //edm4hep::TrackState tsAtFirstHit;
+ //errorCode = marlin_trk.getTrackState(firstHit, tsAtFirstHit, chi2_D, ndf);
+ //debug() << "CombineTracks: TrackState at first hit " << tsAtFirstHit << " chi2 = " << chi2_D << " ndf = " << ndf << endmsg;
+
edm4hep::Vector3d point(0.,0.,0.); // nominal IP
- edm4hep::TrackState trkState ;
- errorCode = marlin_trk.propagate(point, trkState, chi2_D, ndf ) ;
+ edm4hep::TrackState trkState;
+ errorCode = marlin_trk.propagate(point, firstHit, trkState, chi2_D, ndf);
+ debug() << "CombineTracks: output TrackState " << trkState << endmsg;
if ( errorCode != IMarlinTrack::success ) {
debug() << "FullLDCTrackingAlg::CombineTracks: propagate to IP fails with error " << errorCode << endmsg;
@@ -3889,7 +3914,7 @@ void FullLDCTrackingAlg::AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, f
covMatrix[14] = ( _initialTrackError_tanL ); //sigma_tanl^2
pre_fit.covMatrix = covMatrix;
- debug() << "AssignOuterHitsToTracks: before createFit" << endmsg;
+ debug() << "AssignOuterHitsToTracks: before createFit " << pre_fit << endmsg;
int error = MarlinTrk::createFit( trkHits, marlin_trk, &pre_fit, _bField, IMarlinTrack::forward/*backward*/, _maxChi2PerHit );
debug() << "AssignOuterHitsToTracks: after createFit" << endmsg;
if ( error != IMarlinTrack::success ) {
@@ -3912,11 +3937,15 @@ void FullLDCTrackingAlg::AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, f
continue ;
}
debug() << "AssignOuterHitsToTracks: before propagate" << endmsg;
+ std::vector<std::pair<edm4hep::TrackerHit, double> > hitsInFit;
+ marlin_trk->getHitsInFit(hitsInFit);
+ auto firstHit = hitsInFit.begin()->first;
+
edm4hep::Vector3d point(0.,0.,0.); // nominal IP
int return_code = 0;
TrackState trkState ;
- return_code = marlin_trk->propagate(point, trkState, chi2_D, ndf ) ;
+ return_code = marlin_trk->propagate(point, firstHit, trkState, chi2_D, ndf);
delete marlin_trk ;
debug() << "AssignOuterHitsToTracks: after delete" << endmsg;
@@ -4011,7 +4040,8 @@ HelixClass * FullLDCTrackingAlg::GetExtrapolationHelix( TrackExtended * track) {
z_ref_max = fabs(z_ref);
ts_at_calo = getTrackStateAt(trk_lcio, edm4hep::TrackState::AtCalorimeter);
- debug() << "FullLDCTrackingAlg::GetExtrapolationHelix set ts_at_calo with ref_z = " << z_ref << endmsg;
+ debug() << "GetExtrapolationHelix set ts_at_calo with ref_z = " << z_ref << endmsg;
+ debug() << "TrackState at calo: " << ts_at_calo << endmsg;
}
}
}
@@ -4401,7 +4431,7 @@ void FullLDCTrackingAlg::AssignSiHitsToTracks(TrackerHitExtendedVec hitVec,
pre_fit.covMatrix = covMatrix;
- int error = MarlinTrk::createFit( trkHits, marlin_trk, &pre_fit, _bField, IMarlinTrack::backward , _maxChi2PerHit );
+ int error = MarlinTrk::createFit( trkHits, marlin_trk, &pre_fit, _bField, !IMarlinTrack::backward , _maxChi2PerHit );
if ( error != IMarlinTrack::success ) {
debug() << "FullLDCTrackingAlg::AssignSiHitsToTracks: creation of fit fails with error " << error << endmsg;
diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp
index 298861983ba592a17bbf80c9435e3ad634cb1fc8..9de3f481170a69256f9059ac0f2eb73e13d6d69d 100644
--- a/Service/GearSvc/src/GearSvc.cpp
+++ b/Service/GearSvc/src/GearSvc.cpp
@@ -327,7 +327,8 @@ StatusCode GearSvc::convertComposite(dd4hep::DetElement& vtx){
l.zHalfSensitive/dd4hep::mm, l.widthSensitive/dd4hep::mm, 0.);
}
- {
+ if (vtxData->layersPlanar.size()>0) {
+ debug() << "find out planar support material" << endmsg;
const dd4hep::rec::ZPlanarData::LayerLayout& l = vtxData->layersPlanar[0] ;
double offset = l.offsetSupport;
dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, offset, 2.*dd4hep::mm);
@@ -337,7 +338,7 @@ StatusCode GearSvc::convertComposite(dd4hep::DetElement& vtx){
}
std::vector<int> ids;
- std::vector<double> zhalfs, rsens, tsens, rsups, tsups, phi0s, rgaps, dphis;
+ std::vector<double> zhalfs, rsens, tsens, rsups, tsups, phi0s, rgaps, dphis, widths;
for (unsigned i=0,n=vtxData->layersBent.size(); i<n; ++i) {
const dd4hep::rec::CylindricalData::LayerLayout& l = vtxData->layersBent[i];
@@ -351,6 +352,7 @@ StatusCode GearSvc::convertComposite(dd4hep::DetElement& vtx){
phi0s.push_back(l.phi0);
rgaps.push_back(l.rgap/dd4hep::mm);
dphis.push_back(l.dphi);
+ widths.push_back(l.width/dd4hep::mm);
}
gearVTX->setIntVals("VTXLayerIds", ids);
@@ -362,11 +364,13 @@ StatusCode GearSvc::convertComposite(dd4hep::DetElement& vtx){
gearVTX->setDoubleVals("VTXLayerPhi0", phi0s);
gearVTX->setDoubleVals("VTXLayerRadialGap", rgaps);
gearVTX->setDoubleVals("VTXLayerDeltaPhi", dphis);
+ gearVTX->setDoubleVals("VTXLayerWidth", widths);
m_gearMgr->setVXDParameters(gearVTX);
- {
- const dd4hep::rec::CylindricalData::LayerLayout& l = vtxData->layersBent[0] ;
+ if (vtxData->layersBent.size()>0) {
+ debug() << "find out bent support material" << endmsg;
+ 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, "VXDBentSupportMaterial");
@@ -374,6 +378,7 @@ StatusCode GearSvc::convertComposite(dd4hep::DetElement& vtx){
}
if (vtxData->rOuterShell>vtxData->rInnerShell) {
+ debug() << "find out shell material" << endmsg;
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");
@@ -522,7 +527,10 @@ StatusCode GearSvc::convertSIT(dd4hep::DetElement& sit){
sitParams->setDoubleVal("strip_length_mm", sitData->lengthStrip*CLHEP::cm);
sitParams->setDoubleVal("strip_pitch_mm", sitData->pitchStrip*CLHEP::cm);
sitParams->setDoubleVal("strip_angle_deg", strip_angle_deg);
+
std::vector<int> n_sensors_per_ladder;
+ std::vector<double> length_sensors;
+ std::vector<double> thickness_flexs;
for( int layer=0; layer < nLayers; layer++){
dd4hep::rec::ZPlanarData::LayerLayout& layout = sitlayers[layer];
@@ -537,13 +545,50 @@ StatusCode GearSvc::convertSIT(dd4hep::DetElement& sit){
double senOffset = layout.offsetSensitive*CLHEP::cm;
double senThickness = layout.thicknessSensitive*CLHEP::cm;
double senHalfLength = layout.zHalfSensitive*CLHEP::cm;
- double senWidth = layout.widthSensitive*CLHEP::cm;
+ //double senWidth = layout.widthSensitive*CLHEP::cm;
+ // FIXME: treat sensitive width as flex thickness as preliminary, since geometry need same width of support and sensitive
+ double flexThickness = layout.widthSensitive*CLHEP::cm;
int nSensorsPerLadder = layout.sensorsPerLadder;
+ double sensorLength = layout.lengthSensor*CLHEP::cm;
double stripAngle = strip_angle_deg*CLHEP::degree;
n_sensors_per_ladder.push_back(nSensorsPerLadder);
- sitParams->addLayer(nLadders, phi0, supRMin, supOffset, supThickness, supHalfLength, supWidth, 0, senRMin, senOffset, senThickness, senHalfLength, senWidth, 0);
+ length_sensors.push_back(sensorLength);
+ thickness_flexs.push_back(flexThickness);
+
+ if (layer==0) {
+ // support
+ {
+ dd4hep::rec::Vector3D a(layout.distanceSupport, layout.offsetSupport, 0);
+ dd4hep::rec::Vector3D b(layout.distanceSupport + layout.thicknessSupport, layout.offsetSupport, 0);
+ gear::SimpleMaterialImpl* supportMaterial = CreateGearMaterial(a, b, "ITKBarrelSupportMaterial");
+ m_gearMgr->registerSimpleMaterial(supportMaterial);
+ }
+ // sensor
+ {
+ gear::SimpleMaterialImpl* flexMaterial = nullptr;
+ if (senRMin > supRMin) {
+ dd4hep::rec::Vector3D a(layout.distanceSensitive + layout.thicknessSensitive, layout.offsetSensitive, 2.*dd4hep::mm);
+ dd4hep::rec::Vector3D b(layout.distanceSensitive + layout.thicknessSensitive + flexThickness*dd4hep::mm, layout.offsetSensitive, 2.*dd4hep::mm);
+ flexMaterial = CreateGearMaterial(a, b, "ITKBarrelFlexMaterial");
+ }
+ else {
+ dd4hep::rec::Vector3D a(layout.distanceSensitive - flexThickness*dd4hep::mm, layout.offsetSensitive, 2.*dd4hep::mm);
+ dd4hep::rec::Vector3D b(layout.distanceSensitive, layout.offsetSensitive, 2.*dd4hep::mm);
+ flexMaterial = CreateGearMaterial(a, b, "ITKBarrelFlexMaterial");
+ }
+ double ratio = 1.0 + flexThickness/flexMaterial->getRadLength()*93.6607/senThickness;
+ debug() << "sensor thickness: " << senThickness << " flex thickness: " << flexThickness << " radL: " << flexMaterial->getRadLength() << " effetive thickness: " << ratio*senThickness << endmsg;
+
+ gear::SimpleMaterialImpl* senMaterial = new gear::SimpleMaterialImpl("ITKBarrelSensorMaterial", 28.09, 14.0, 2330.*ratio, 93.6607/ratio, 0);
+ m_gearMgr->registerSimpleMaterial(flexMaterial);
+ m_gearMgr->registerSimpleMaterial(senMaterial);
+ }
+ }
+ sitParams->addLayer(nLadders, phi0, supRMin, supOffset, supThickness, supHalfLength, supWidth, 0, senRMin, senOffset, senThickness, senHalfLength, supWidth, 0);
}
sitParams->setIntVals("n_sensors_per_ladder",n_sensors_per_ladder);
+ sitParams->setDoubleVals("length_sensors", length_sensors);
+ sitParams->setDoubleVals("thickness_flexs", thickness_flexs);
m_gearMgr->setSITParameters( sitParams ) ;
return StatusCode::SUCCESS;
@@ -589,7 +634,28 @@ 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 );
+ // InnerWall
+ {
+ dd4hep::rec::Vector3D a(tpcData->rMin, 0, tpcData->driftLength/2.0);
+ dd4hep::rec::Vector3D b(tpcData->rMin+tpcData->innerWallThickness, 0, tpcData->driftLength/2.0);
+ gear::SimpleMaterialImpl* TPCInnerWallMaterial = CreateGearMaterial(a, b, "TPCInnerWallMaterial");
+ m_gearMgr->registerSimpleMaterial(TPCInnerWallMaterial);
+ }
double x = (supportData->sections[0].rInner + supportData->sections[0].rOuter)/2.;
+ // Gas
+ {
+ dd4hep::rec::Vector3D a(x, 0, tpcData->driftLength/2.0);
+ dd4hep::rec::Vector3D b(x+1*dd4hep::mm, 0, tpcData->driftLength/2.0);
+ gear::SimpleMaterialImpl* TPCGasMaterial = CreateGearMaterial(a, b, "TPCGasMaterial");
+ m_gearMgr->registerSimpleMaterial(TPCGasMaterial);
+ }
+ // OuterWall
+ {
+ dd4hep::rec::Vector3D a(tpcData->rMax-tpcData->outerWallThickness, 0, tpcData->driftLength/2.0);
+ dd4hep::rec::Vector3D b(tpcData->rMax, 0, tpcData->driftLength/2.0);
+ gear::SimpleMaterialImpl* TPCOuterWallMaterial = CreateGearMaterial(a, b, "TPCOuterWallMaterial");
+ m_gearMgr->registerSimpleMaterial(TPCOuterWallMaterial);
+ }
// Readout
{
dd4hep::rec::Vector3D a(x, 0, supportData->sections[0].zPos);
@@ -790,6 +856,16 @@ StatusCode GearSvc::convertSET(dd4hep::DetElement& set){
int nSensorsPerLadder = layout.sensorsPerLadder;
double stripAngle = strip_angle_deg*CLHEP::degree;
n_sensors_per_ladder.push_back(nSensorsPerLadder);
+
+ if (layer==0) {
+ // support
+ {
+ dd4hep::rec::Vector3D a(layout.distanceSupport, layout.offsetSupport, 0);
+ dd4hep::rec::Vector3D b(layout.distanceSupport + layout.thicknessSupport, layout.offsetSupport, 0);
+ gear::SimpleMaterialImpl* supportMaterial = CreateGearMaterial(a, b, "OTKBarrelSupportMaterial");
+ m_gearMgr->registerSimpleMaterial(supportMaterial);
+ }
+ }
setParams->addLayer(nLadders, phi0, supRMin, supOffset, supThickness, supHalfLength, supWidth, 0, senRMin, senOffset, senThickness, senHalfLength, senWidth, 0);
}
setParams->setIntVals("n_sensors_per_ladder", n_sensors_per_ladder);
@@ -891,7 +967,9 @@ gear::SimpleMaterialImpl* GearSvc::CreateGearMaterial(const dd4hep::rec::Vector3
// TODO: move to GeomSvc
dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() );
- const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween(a, b);
+ // 0.0005 mm make sure don't loss material because of precision
+ dd4hep::rec::Vector3D safe_env = (a.z() == b.z()) ? dd4hep::rec::Vector3D(0.0005*dd4hep::mm, 0, 0) : dd4hep::rec::Vector3D(0, 0, 0.0005*dd4hep::mm);
+ const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween(a-safe_env, b+safe_env);
dd4hep::rec::MaterialData mat = (materials.size() > 1) ? matMgr.createAveragedMaterial(materials) : materials[0].first;
debug() << " ####### found materials between points : " << a << " and " << b << " ######" << endmsg;
diff --git a/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h b/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h
index 70303ae34a3c9547e2cfbe359e7cc651aace6287..c008cbc914bd8592c2d36b0cee674d306700fc6c 100644
--- a/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h
+++ b/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h
@@ -74,6 +74,7 @@ namespace MarlinTrk{
IMarlinTrack* marlinTrk,
edm4hep::MutableTrack* track,
std::vector<edm4hep::TrackerHit>& hit_list,
+ bool fit_backwards,
edm4hep::TrackState* atLastHit=0,
edm4hep::TrackState* atCaloFace=0);
@@ -82,7 +83,8 @@ namespace MarlinTrk{
/** Set the subdetector hit numbers for the TrackImpl */
void addHitNumbersToTrack(edm4hep::MutableTrack* track, std::vector<std::pair<edm4hep::TrackerHit , double> >& hit_list, bool hits_in_fit, UTIL::BitField64& cellID_encoder);
-
+
+ int createTrackStateAtCaloFace( IMarlinTrack* marlinTrk, edm4hep::TrackState* track, edm4hep::TrackerHit trkhit, bool tanL_is_positive );
}
#endif
diff --git a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
index 2a4c9d8a8bffc5ca91f03c57d91fc6eb9c3451df..0a4f6bf6c0a89a44b5c92708456155e863ad0bd9 100644
--- a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
+++ b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc
@@ -103,7 +103,7 @@ namespace MarlinTrk {
- int createTrackStateAtCaloFace( IMarlinTrack* marlinTrk, edm4hep::TrackState* track, edm4hep::TrackerHit trkhit, bool tanL_is_positive );
+ //int createTrackStateAtCaloFace( IMarlinTrack* marlinTrk, edm4hep::TrackState* track, edm4hep::TrackerHit trkhit, bool tanL_is_positive );
int createFinalisedLCIOTrack( IMarlinTrack* marlinTrk, std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::MutableTrack* track, bool fit_backwards, const decltype(edm4hep::TrackState::covMatrix)& initial_cov_for_prefit, float bfield_z, double maxChi2Increment){
@@ -166,7 +166,7 @@ namespace MarlinTrk {
return fit_status;
}
- int error = finaliseLCIOTrack(marlinTrk, track, hit_list);
+ int error = finaliseLCIOTrack(marlinTrk, track, hit_list, fit_backwards);
//std::cout << "finaliseLCIOTrack. status = " << error << std::endl;
return error;
}
@@ -212,16 +212,7 @@ namespace MarlinTrk {
///////////////////////////////////////////////////////
// add hits to IMarlinTrk
///////////////////////////////////////////////////////
- std::vector<edm4hep::TrackerHit> hit_list_copy;
- if (fit_backwards) {
- std::reverse_copy(hit_list.begin(), hit_list.end(), std::back_inserter(hit_list_copy));
- }
- else {
- hit_list_copy.reserve(hit_list.size());
- hit_list_copy.insert(hit_list_copy.end(), hit_list.begin(), hit_list.end());
- }
- hit_list_copy.swap(hit_list);
std::vector<edm4hep::TrackerHit>::iterator it = hit_list.begin();
// start by trying to add the hits to the track we want to finally use.
@@ -260,9 +251,10 @@ namespace MarlinTrk {
if (isSuccessful) {
added_hits.push_back(trkHit);
+ //std::cout << "DEBUG<<<<<MarlinTrkUtils::createFit Hit " << it - hit_list.begin() << " Added " << trkHit.getPosition() << std::endl;
}
else{
- //std::cout << "DEBUG<<<<<MarlinTrkUtils::createFit Hit " << it - hit_list.begin() << " Dropped " << std::endl;
+ //std::cout << "DEBUG<<<<<MarlinTrkUtils::createFit Hit " << it - hit_list.begin() << " Dropped " << trkHit.getPosition() << std::endl;
}
}
@@ -272,7 +264,6 @@ namespace MarlinTrk {
return IMarlinTrack::bad_intputs;
}
- hit_list_copy.swap(hit_list);
///////////////////////////////////////////////////////
// set the initial track parameters
///////////////////////////////////////////////////////
@@ -344,16 +335,18 @@ namespace MarlinTrk {
const edm4hep::Vector3d& x1 = twoD_hits[0].getPosition();
const edm4hep::Vector3d& x2 = twoD_hits[ twoD_hits.size()/2 ].getPosition();
const edm4hep::Vector3d& x3 = twoD_hits.back().getPosition();
-
+
HelixTrack helixTrack( x1, x2, x3, bfield_z, HelixTrack::forwards );
- if ( fit_backwards == IMarlinTrack::backward ) {
- pre_fit->location = MarlinTrk::Location::AtLastHit;
- helixTrack.moveRefPoint(hit_list.back().getPosition()[0], hit_list.back().getPosition()[1], hit_list.back().getPosition()[2]);
- } else {
- pre_fit->location = MarlinTrk::Location::AtFirstHit;
- helixTrack.moveRefPoint(hit_list.front().getPosition()[0], hit_list.front().getPosition()[1], hit_list.front().getPosition()[2]);
- }
+ helixTrack.moveRefPoint(0.0, 0.0, 0.0);
+
+ //if ( fit_backwards == IMarlinTrack::backward ) {
+ //pre_fit->location = MarlinTrk::Location::AtLastHit;
+ //helixTrack.moveRefPoint(hit_list.back().getPosition()[0], hit_list.back().getPosition()[1], hit_list.back().getPosition()[2]);
+ //} else {
+ //pre_fit->location = MarlinTrk::Location::AtFirstHit;
+ //helixTrack.moveRefPoint(hit_list.front().getPosition()[0], hit_list.front().getPosition()[1], hit_list.front().getPosition()[2]);
+ //}
const float referencePoint[3] = { helixTrack.getRefPointX() , helixTrack.getRefPointY() , helixTrack.getRefPointZ() };
@@ -369,7 +362,7 @@ namespace MarlinTrk {
}
- int finaliseLCIOTrack( IMarlinTrack* marlintrk, edm4hep::MutableTrack* track, std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::TrackState* atLastHit, edm4hep::TrackState* atCaloFace){
+ int finaliseLCIOTrack( IMarlinTrack* marlintrk, edm4hep::MutableTrack* track, std::vector<edm4hep::TrackerHit>& hit_list, bool fit_backwards, edm4hep::TrackState* atLastHit, edm4hep::TrackState* atCaloFace){
///////////////////////////////////////////////////////
// check inputs
@@ -494,7 +487,7 @@ namespace MarlinTrk {
///////////////////////////////////////////////////////
edm4hep::TrackState* trkStateAtFirstHit = new edm4hep::TrackState() ;
- edm4hep::TrackerHit firstHit = hits_in_fit.back().first;
+ edm4hep::TrackerHit firstHit = ( fit_backwards == IMarlinTrack::backward ? hits_in_fit.back().first : hits_in_fit.front().first ) ;
///////////////////////////////////////////////////////
// last hit
@@ -502,7 +495,7 @@ namespace MarlinTrk {
edm4hep::TrackState* trkStateAtLastHit = new edm4hep::TrackState() ;
- edm4hep::TrackerHit lastHit = hits_in_fit.front().first;
+ edm4hep::TrackerHit lastHit = ( fit_backwards == IMarlinTrack::backward ? hits_in_fit.front().first : hits_in_fit.back().first ) ;
#if PODIO_BUILD_VERSION < PODIO_VERSION(0, 17, 4)
edm4hep::TrackerHit last_constrained_hit(0);// = 0 ;
@@ -532,9 +525,18 @@ namespace MarlinTrk {
///////////////////////////////////////////////////////
// make sure that the track state can be propagated to the IP
///////////////////////////////////////////////////////
-
- return_error = marlintrk->propagate(point, firstHit, *trkStateIP, chi2, ndf ) ;
-
+
+ //FIXME: AidaTT not used now, if to use AidaTT, these code need to update and apply
+ //MarlinTrk::IMarlinTrkSystem* trksystem = MarlinTrk::Factory::getCurrentMarlinTrkSystem() ;
+ //bool usingAidaTT = ( trksystem->name() == "AidaTT" ) ;
+ bool usingAidaTT = false;
+ if (fit_backwards == IMarlinTrack::backward || usingAidaTT ) {
+ return_error = marlintrk->propagate(point, firstHit, *trkStateIP, chi2, ndf ) ;
+ }
+ else {
+ //FIXME: if forward, better by add the last inner hits with a Kalman step from the last_constrained hit
+ return_error = marlintrk->propagate(point, firstHit, *trkStateIP, chi2, ndf ) ;
+ }
if ( return_error != IMarlinTrack::success ) {
//streamlog_out(DEBUG4) << "MarlinTrk::finaliseLCIOTrack: return_code for propagation = " << return_error << " NDF = " << ndf << std::endl;
delete trkStateIP;
@@ -613,6 +615,7 @@ namespace MarlinTrk {
if ( return_error == 0 ) {
//std::cout << "fucdout referencePoint " << trkStateCalo.referencePoint << std::endl;
trkStateCalo.location = MarlinTrk::Location::AtCalorimeter;
+ // std::cout << "createTrackStateAtCaloFace OK: " << trkStateCalo << std::endl;
track->addToTrackStates(trkStateCalo);
} else {
#ifdef DEBUG
@@ -663,6 +666,8 @@ namespace MarlinTrk {
}
int return_error = 0;
+ int return_error_barrel = 0;
+ int return_error_endcap = 0;
double chi2 = -DBL_MAX;
int ndf = 0;
@@ -670,28 +675,66 @@ namespace MarlinTrk {
UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string );
encoder.reset() ; // reset to 0
- encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::ECAL ;
- encoder[lcio::ILDCellID0::side] = lcio::ILDDetID::barrel;
- encoder[lcio::ILDCellID0::layer] = 0 ;
+ // ================== need to get the correct ID(s) for the calorimeter face ============================
+ unsigned ecal_barrel_face_ID = lcio::ILDDetID::ECAL ;
+ unsigned ecal_endcap_face_ID = lcio::ILDDetID::ECAL_ENDCAP ;
+
+ //=========================================================================================================
+
+ encoder[lcio::LCTrackerCellID::subdet()] = ecal_barrel_face_ID ;
+ encoder[lcio::LCTrackerCellID::side()] = lcio::ILDDetID::barrel;
+ encoder[lcio::LCTrackerCellID::layer()] = 0 ;
int detElementID = 0;
-
- return_error = marlintrk->propagateToLayer(encoder.lowWord(), trkhit, *trkStateCalo, chi2, ndf, detElementID, IMarlinTrack::modeForward ) ;
-
- if (return_error == IMarlinTrack::no_intersection ) { // try forward or backward
- if (tanL_is_positive) {
- encoder[lcio::ILDCellID0::side] = lcio::ILDDetID::fwd;
+ edm4hep::TrackState tsBarrel;
+ edm4hep::TrackState tsEndcap;
+ // propagate to the barrel layer
+ return_error_barrel = marlintrk->propagateToLayer(encoder.lowWord(), trkhit, tsBarrel, chi2, ndf, detElementID, IMarlinTrack::modeForward ) ;
+ // propagate to the endcap layer
+ encoder[lcio::LCTrackerCellID::subdet()] = ecal_endcap_face_ID ;
+ if (tanL_is_positive) {
+ encoder[lcio::LCTrackerCellID::side()] = lcio::ILDDetID::fwd;
+ }
+ else{
+ encoder[lcio::LCTrackerCellID::side()] = lcio::ILDDetID::bwd;
+ }
+ return_error_endcap = marlintrk->propagateToLayer(encoder.lowWord(), trkhit, tsEndcap, chi2, ndf, detElementID, IMarlinTrack::modeForward ) ;
+
+ // check which is the right intersection / closer to the trkhit
+ if ( return_error_barrel == IMarlinTrack::no_intersection ){
+ //if barrel fails just return ts at the Endcap if exists
+ return_error = return_error_endcap;
+ *trkStateCalo = tsEndcap;
+ }
+ else if( return_error_endcap == IMarlinTrack::no_intersection ){
+ //if barrel succeeded and endcap fails return ts at the barrel
+ return_error = return_error_barrel;
+ *trkStateCalo = tsBarrel;
+ }
+ else{
+ //this means both barrel and endcap have intersections. Return closest to the tracker hit
+ auto hitPos = trkhit.getPosition();
+ auto barrelPos = tsBarrel.referencePoint;
+ auto endcapPos = tsEndcap.referencePoint;
+ double dToBarrel = sqrt((hitPos.x-barrelPos.x)*(hitPos.x-barrelPos.x)+(hitPos.y-barrelPos.y)*(hitPos.y-barrelPos.y)+(hitPos.z-barrelPos.z)*(hitPos.z-barrelPos.z));
+ double dToendcap = sqrt((hitPos.x-endcapPos.x)*(hitPos.x-endcapPos.x)+(hitPos.y-endcapPos.y)*(hitPos.y-endcapPos.y)+(hitPos.z-endcapPos.z)*(hitPos.z-endcapPos.z));
+ if ( dToBarrel < dToendcap ){
+ return_error = return_error_barrel;
+ *trkStateCalo = tsBarrel;
}
else{
- encoder[lcio::ILDCellID0::side] = lcio::ILDDetID::bwd;
+ return_error = return_error_endcap;
+ *trkStateCalo = tsEndcap;
}
- return_error = marlintrk->propagateToLayer(encoder.lowWord(), trkhit, *trkStateCalo, chi2, ndf, detElementID, IMarlinTrack::modeForward ) ;
}
if (return_error !=IMarlinTrack::success ) {
+ //std::cout << "createTrackStateAtCaloFace : could not get TrackState at Calo Face: return_error = " << return_error_barrel << " " << return_error_endcap << std::endl ;
//streamlog_out( DEBUG5 ) << " >>>>>>>>>>> createTrackStateAtCaloFace : could not get TrackState at Calo Face: return_error = " << return_error << std::endl ;
}
-
+ else {
+ //std::cout << "createTrackStateAtCaloFace : " << (*trkStateCalo) << std::endl;
+ }
return return_error;
diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
index 7fcea7993007df90b8634f9f85554f6da2748fe3..3850dc2f1d6a3be69c172eaeda8413c205cb2e45 100644
--- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
+++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
@@ -64,7 +64,7 @@ private:
Gaudi::Property<std::vector<std::string>> m_trackerColNames{this,
"TrackerCollections",
{"VXD", "FTD", "SIT", "TPC", "TPCLowPt", "TPCSpacePoint", "SET",
- "OTKBarrel", "OTKEndcap", "COIL", "MuonBarrel", "MuonEndcap"},
+ "ITKBarrel", "OTKBarrel", "OTKEndcap", "COIL", "MuonBarrel", "MuonEndcap"},
"Names of the Tracker collections (without suffix Collection)"};
Gaudi::Property<std::vector<std::string>> m_calorimeterColNames{this,
"CalorimeterCollections",
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
index 94329ee4fff90dd10cca21d83d2711509d21add1..e75bd1ded90a7130f49b6198af80fd917e34a880 100644
--- a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
+++ b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
@@ -26,7 +26,8 @@ G4bool GenericTrackerSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHis
G4TouchableHandle touchPost = step->GetPostStepPoint()->GetTouchableHandle();
G4TouchableHandle touchPre = step->GetPreStepPoint()->GetTouchableHandle();
dd4hep::sim::Geant4StepHandler h(step);
- if (fabs(h.trackDef()->GetPDGCharge()) < 0.01) return true;
+ //if (fabs(h.trackDef()->GetPDGCharge()) < 0.01) return true;
+ if (h.deposit() <= 0 && !Geant4Hit::isGeantino(step->GetTrack())) return true;
dd4hep::Position prePos = h.prePos();
dd4hep::Position postPos = h.postPos();
diff --git a/Utilities/KalDet/include/kaldet/CEPCArcDiscMeasLayer.h b/Utilities/KalDet/include/kaldet/CEPCArcDiscMeasLayer.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed8ad0fa4658f871dc5b33cf7b0648ffbc4456a1
--- /dev/null
+++ b/Utilities/KalDet/include/kaldet/CEPCArcDiscMeasLayer.h
@@ -0,0 +1,115 @@
+#ifndef __CEPCARCDISCMEASLAYER__
+#define __CEPCARCDISCMEASLAYER__
+
+/** CEPCArcDiscMeasLayer: User defined KalTest Disc measurement layer class used with ILDPLanarTrackHit.
+ *
+ * @author S.Aplin DESY
+ */
+
+#include "TVector3.h"
+
+#include "kaltest/TKalMatrix.h"
+#include "kaltest/TPlane.h"
+#include "kaltest/KalTrackDim.h"
+#include "ILDVMeasLayer.h"
+
+#include "TMath.h"
+#include <sstream>
+
+class TVTrackHit;
+
+
+class CEPCArcDiscMeasLayer : public ILDVMeasLayer, public TPlane {
+
+public:
+
+ /** Constructor Taking inner and outer materials, center and normal to the plane, B-Field, sorting policy, min and max r, whether the layer is sensitive, Cell ID, and an optional name */
+
+ CEPCArcDiscMeasLayer(TMaterial &min,
+ TMaterial &mout,
+ const TVector3 ¢er,
+ const TVector3 &normal,
+ double Bz,
+ double SortingPolicy,
+ double rMin,
+ double rMax,
+ Bool_t is_active,
+ Int_t CellID = -1,
+ const Char_t *name = "CEPCArcDiscMeasL")
+ : ILDVMeasLayer(min, mout, Bz, is_active, CellID, name),
+ TPlane(center, normal),
+ _sortingPolicy(SortingPolicy), _rMin(rMin), _rMax(rMax)
+ { /* no op */ }
+
+ CEPCArcDiscMeasLayer(TMaterial &min,
+ TMaterial &mout,
+ const TVector3 ¢er,
+ const TVector3 &normal,
+ double Bz,
+ double SortingPolicy,
+ double rMin,
+ double rMax,
+ const std::vector<int>& CellIDs,
+ Bool_t is_active,
+ const Char_t *name = "CEPCArcDiscMeasL")
+ : ILDVMeasLayer(min, mout, Bz, CellIDs, is_active, name),
+ TPlane(center, normal),
+ _sortingPolicy(SortingPolicy), _rMin(rMin), _rMax(rMax)
+ { /* no op */ }
+
+ // Parrent's pure virtuals that must be implemented
+
+ /** Global to Local coordinates */
+ virtual TKalMatrix XvToMv (const TVTrackHit &ht,
+ const TVector3 &xv) const
+ { return this->XvToMv(xv); }
+
+ /** Global to Local coordinates */
+ virtual TKalMatrix XvToMv (const TVector3 &xv) const;
+
+ /** Local to Global coordinates */
+ virtual TVector3 HitToXv (const TVTrackHit &ht) const;
+
+ /** Calculate Projector Matrix */
+ virtual void CalcDhDa (const TVTrackHit &ht,
+ const TVector3 &xv,
+ const TKalMatrix &dxphiada,
+ TKalMatrix &H) const;
+
+ virtual Int_t CalcXingPointWith(const TVTrack &hel,
+ TVector3 &xx,
+ Double_t &phi,
+ Int_t mode,
+ Double_t eps) const;
+
+
+ /** Convert LCIO Tracker Hit to an ILDPLanarTrackHit */
+ virtual ILDVTrackHit* ConvertLCIOTrkHit(edm4hep::TrackerHit trkhit) const ;
+
+ /** Check if global point is on surface */
+ inline virtual Bool_t IsOnSurface (const TVector3 &xx) const;
+
+ /** Get sorting policy for this plane */
+ double GetSortingPolicy() const { return _sortingPolicy; }
+
+ /** Get the intersection and the CellID, needed for multilayers */
+ virtual int getIntersectionAndCellID(const TVTrack &hel,
+ TVector3 &xx,
+ Double_t &phi,
+ Int_t &CellID,
+ Int_t mode,
+ Double_t eps = 1.e-8) const {
+
+ CellID = this->getCellIDs()[0]; // not multilayer
+ return this->CalcXingPointWith(hel,xx,phi,0,eps);
+
+ }
+
+private:
+ double _sortingPolicy;
+ double _rMin;
+ double _rMax;
+
+};
+
+#endif
diff --git a/Utilities/KalDet/include/kaldet/CEPCCylinderMeasLayer.h b/Utilities/KalDet/include/kaldet/CEPCCylinderMeasLayer.h
new file mode 100644
index 0000000000000000000000000000000000000000..619ac756ff12b9fff9359d83099176add2625f79
--- /dev/null
+++ b/Utilities/KalDet/include/kaldet/CEPCCylinderMeasLayer.h
@@ -0,0 +1,57 @@
+#ifndef CEPCCYLINDERMEASLAYER_H
+#define CEPCCYLINDERMEASLAYER_H
+
+/** CEPCCylinderMeasLayer: User defined KalTest measurement layer class
+ *
+ * @author
+ */
+
+#include "ILDCylinderMeasLayer.h"
+#include <iostream>
+#include <cmath>
+/* #include "streamlog/streamlog.h" */
+
+
+class CEPCCylinderMeasLayer : public ILDCylinderMeasLayer/*, public TCylinder */{
+
+ public:
+
+ /** Constructor Taking inner and outer materials, radius and half length, B-Field, whether the layer is sensitive, Cell ID, and an optional name */
+ CEPCCylinderMeasLayer(TMaterial &min,
+ TMaterial &mout,
+ Double_t r0,
+ Double_t lhalf,
+ Double_t phi0,
+ Double_t width,
+ Double_t x0,
+ Double_t y0,
+ Double_t z0,
+ Double_t Bz,
+ Bool_t is_active,
+ Int_t CellID = -1,
+ const Char_t *name = "CEPCCylinderMeasL");
+ // : ILDCylinderMeasLayer(min, mout, r0, lhalf, x0, y0, z0, Bz, is_active, CellID, name),
+ //_phi0(phi0), _dphi(dphi)
+ //{ /* no op */ }
+
+ virtual Bool_t IsOnSurface(const TVector3 &xx) const;// {
+
+ // bool z = (xx.Z() >= GetZmin() && xx.Z() <= GetZmax());
+ // bool r = std::fabs( (xx-this->GetXc()).Perp() - this->GetR() ) < 1.e-3; // for very short, very stiff tracks this can be poorly defined, so we relax this here a bit to 1 micron
+ // Double_t phi = (xx-this->GetXc()).Phi();
+ //dphi = phi - phi0;
+ //if (dphi>M_PI) dphi -= M_PI*2;
+ //else if (dphi<-M_PI) dphi += M_PI*2;
+// streamlog_out(DEBUG0) << "CEPCCylinderMeasLayer IsOnSurface for " << this->TVMeasLayer::GetName() << " R = " << this->GetR() << " GetZmin() = " << GetZmin() << " GetZmax() = " << GetZmax()
+// << " dr = " << std::fabs( (xx-this->GetXc()).Perp() - this->GetR() ) << " r = " << r << " z = " << z
+// << std::endl;
+
+ //return r && z && (fabs(dphi)<_dphi);
+ //}
+
+private:
+
+ Double_t _phi0;
+ Double_t _dphi;
+};
+#endif
diff --git a/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h b/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h
index 7184f78c5c8b25ab1b2373fa3be15cf7f36a5b38..9da9598df62f7a3dac0bb5c5d3d24c0f919c3139 100644
--- a/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h
+++ b/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h
@@ -52,6 +52,7 @@ private:
struct STT_Layer {
int id;
double length;
+ double width;
double senRMin;
double senThickness;
double supRMin;
diff --git a/Utilities/KalDet/include/kaldet/ILDCylinderMeasLayer.h b/Utilities/KalDet/include/kaldet/ILDCylinderMeasLayer.h
index 00914f6580882614e07115bfe1e1fb3ad13f33ca..08a3c82c0b80ef34550a6045f417d06fe47cb89e 100644
--- a/Utilities/KalDet/include/kaldet/ILDCylinderMeasLayer.h
+++ b/Utilities/KalDet/include/kaldet/ILDCylinderMeasLayer.h
@@ -34,7 +34,7 @@ public:
{ /* no op */ }
- Bool_t IsOnSurface(const TVector3 &xx) const {
+ virtual Bool_t IsOnSurface(const TVector3 &xx) const {
bool z = (xx.Z() >= GetZmin() && xx.Z() <= GetZmax());
bool r = std::fabs( (xx-this->GetXc()).Perp() - this->GetR() ) < 1.e-3; // for very short, very stiff tracks this can be poorly defined, so we relax this here a bit to 1 micron
diff --git a/Utilities/KalDet/src/ild/common/CEPCArcDiscMeasLayer.cc b/Utilities/KalDet/src/ild/common/CEPCArcDiscMeasLayer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..edf0c61a94aa0042dd2b582fb0af6be3224afca1
--- /dev/null
+++ b/Utilities/KalDet/src/ild/common/CEPCArcDiscMeasLayer.cc
@@ -0,0 +1,301 @@
+
+#include <iostream>
+
+#include "kaldet/CEPCArcDiscMeasLayer.h"
+#include "kaldet/ILDPlanarHit.h"
+
+#include "kaltest/TVTrack.h"
+#include "TVector3.h"
+#include "TMath.h"
+#include "TRotMatrix.h"
+#include "TBRIK.h"
+#include "TNode.h"
+#include "TString.h"
+
+#include <edm4hep/TrackerHit.h>
+
+#include "gearimpl/Vector3D.h"
+
+#include "DetIdentifier/CEPCConf.h"
+#include <bitset>
+// #include "streamlog/streamlog.h"
+
+
+TKalMatrix CEPCArcDiscMeasLayer::XvToMv(const TVector3 &xv) const
+{
+
+ // Calculate measurement vector (hit coordinates) from global coordinates:
+
+ TKalMatrix mv(kMdim,1);
+
+ mv(0,0) = xv.X() ;
+
+
+ mv(1,0) = xv.Y() ;
+ return mv;
+
+}
+
+
+TVector3 CEPCArcDiscMeasLayer::HitToXv(const TVTrackHit &vht) const
+{
+ // const ILDPlanarHit &mv = dynamic_cast<const ILDPlanarHit &>(vht);
+
+ double x = vht(0,0) ;
+ double y = vht(1,0) ;
+
+ double z = GetXc().Z() ;
+
+ return TVector3(x,y,z);
+}
+
+void CEPCArcDiscMeasLayer::CalcDhDa(const TVTrackHit &vht,
+ const TVector3 &xxv,
+ const TKalMatrix &dxphiada,
+ TKalMatrix &H) const
+{
+ // Calculate
+ // H = (@h/@a) = (@phi/@a, @z/@a)^t
+ // where
+ // h(a) = (phi, z)^t: expected meas vector
+ // a = (drho, phi0, kappa, dz, tanl, t0)
+ //
+
+ Int_t sdim = H.GetNcols();
+ Int_t hdim = TMath::Max(5,sdim-1);
+
+ // Set H = (@h/@a) = (@d/@a, @z/@a)^t
+
+ for (Int_t i=0; i<hdim; i++) {
+
+ H(0,i) = dxphiada(0,i);
+ H(1,i) = dxphiada(1,i) ;
+
+ }
+ if (sdim == 6) {
+ H(0,sdim-1) = 0.;
+ H(1,sdim-1) = 0.;
+ }
+
+}
+
+Int_t CEPCArcDiscMeasLayer::CalcXingPointWith(const TVTrack &hel,
+ TVector3 &xx,
+ Double_t &phi,
+ Int_t mode,
+ Double_t eps) const{
+
+ phi = 0.0;
+
+ xx.SetX(0.0);
+ xx.SetY(0.0);
+ xx.SetZ(0.0);
+
+
+ // check that direction has one of the correct values
+ if( !( mode == 0 || mode == 1 || mode == -1) ) return -1 ;
+
+ // get helix parameters
+ Double_t dr = hel.GetDrho();
+ Double_t phi0 = hel.GetPhi0(); //
+ Double_t kappa = hel.GetKappa();
+ Double_t rho = hel.GetRho();
+ Double_t omega = 1.0 / rho;
+ Double_t z0 = hel.GetDz();
+ Double_t tanl = hel.GetTanLambda();
+
+ TVector3 ref_point = hel.GetPivot();
+
+
+ //
+ // Check if charge is nonzero.
+ //
+
+ Int_t chg = (Int_t)TMath::Sign(1.1,kappa);
+ if (!chg) {
+ // streamlog_out(ERROR) << ">>>> Error >>>> CEPCArcDiscMeasLayer::CalcXingPointWith" << std::endl
+ // << " Kappa = 0 is invalid for a helix " << std::endl;
+ return -1;
+ }
+
+ const double sin_phi0 = sin(phi0);
+ const double cos_phi0 = cos(phi0);
+
+ const double x_pca = ref_point.x() + dr * cos_phi0 ;
+ const double y_pca = ref_point.y() + dr * sin_phi0 ;
+ const double z_pca = ref_point.z() + z0 ;
+
+ const double z = this->GetXc().Z() ;
+ // get path length to crossing point
+
+ const double s = ( z - z_pca ) / tanl ;
+
+// streamlog_out(DEBUG0) << "CEPCArcDiscMeasLayer::CalcXingPointWith "
+// << " ref_point.z() = " << ref_point.z()
+// << " z = " << z
+// << " z0 = " << z0
+// << " z_pca = " << z_pca
+// << " tanl = " << tanl
+// << " z - z_pca = " << z - z_pca
+// << std::endl;
+
+// TVector3 xx_n;
+// int cuts = TVSurface::CalcXingPointWith(hel, xx_n, phi, 0, eps);
+// streamlog_out(DEBUG0) << "CEPCArcDiscMeasLayer::CalcXingPointWith from Newton: cuts = " << cuts << " x = " << xx_n.x() << " y = "<< xx_n.y() << " z = " << xx_n.z() << " r = " << xx_n.Perp() << " phi = " << xx_n.Phi() << " dphi = " << phi << std::endl;
+
+
+ phi = -omega * s;
+
+ const double delta_phi_half = -phi/2.0 ;
+
+
+ double x;
+ double y;
+
+ if( fabs(s) > FLT_MIN ){ // protect against starting on the plane
+
+ x = x_pca - s * ( sin(delta_phi_half) / delta_phi_half ) * sin( phi0 - delta_phi_half ) ;
+
+ y = y_pca + s * ( sin(delta_phi_half) / delta_phi_half ) * cos( phi0 - delta_phi_half ) ;
+
+ }
+ else{
+ // streamlog_out(DEBUG0) << "CEPCArcDiscMeasLayer::CalcXingPointWith Using PCA values " << std::endl;
+ x = x_pca;
+ y = y_pca;
+ phi = 0;
+ }
+
+
+ // check if intersection with plane is within boundaries
+
+ xx.SetXYZ(x, y, z);
+
+
+ // streamlog_out(DEBUG0) << "CEPCArcDiscMeasLayer::CalcXingPointWith : cuts = " << (IsOnSurface(xx) ? 1 : 0) << " x = " << xx.x() << " y = "<< xx.y() << " z = " << xx.z() << " r = " << xx.Perp() << " phi = " << xx.Phi() << " dphi = " << phi << " s = " << s << " " << this->TVMeasLayer::GetName() << std::endl;
+
+ if( mode!=0 && fabs(phi)>1.e-10 ){ // (+1,-1) = (fwd,bwd)
+ if( chg*phi*mode > 0){
+ return 0;
+ }
+ }
+
+ return (IsOnSurface(xx) ? 1 : 0);
+
+}
+
+
+Bool_t CEPCArcDiscMeasLayer::IsOnSurface(const TVector3 &xx) const
+{
+
+ bool onSurface = false ;
+
+ TKalMatrix mv = XvToMv(xx);
+
+ // check whether the hit lies in the same plane as the surface
+ if (TMath::Abs((xx.X()-GetXc().X())*GetNormal().X() + (xx.Y()-GetXc().Y())*GetNormal().Y() + (xx.Z()-GetXc().Z())*GetNormal().Z()) < 1e-4) {
+ // check whether the hit lies within the boundary of the surface
+
+ double r2 = mv(0,0) * mv(0,0) + mv(1,0) * mv(1,0) ;
+
+ if (r2 <= _rMax*_rMax && r2 >= _rMin*_rMin) {
+ onSurface = true ;
+ }
+ else {
+ std::cout << "r2: " << r2 << " r2min: " << _rMin*_rMin << " r2max: " << _rMax*_rMax << std::endl;
+ }
+ }
+ else {
+ std::cout << "Xc: " << GetXc().X() << " " << GetXc().Y() << " " << GetXc().Z() << " Normal: " << GetNormal().X() << " " << GetNormal().Y() << " " << GetNormal().Z() << std::endl;
+ }
+
+ return onSurface;
+
+}
+
+
+ILDVTrackHit* CEPCArcDiscMeasLayer::ConvertLCIOTrkHit(edm4hep::TrackerHit trkhit) const {
+
+ //edm4hep::TrackerHitPlane* plane_hit = dynamic_cast<EVENT::TrackerHitPlane*>( trkhit ) ;
+ //edm4hep::TrackerHitPlane* plane_hit = trkhit;
+ std::cout << "CEPCArcDiscMeasLayer::ConvertLCIOTrkHit type = " << trkhit.getType() << std::endl;
+ std::bitset<32> type(trkhit.getType());
+ //if (!type[CEPCConf::TrkHitTypeBit::PLANAR]) return NULL;
+
+ //edm4hep::TrackerHit plane_hit = trkhit;
+ //if( plane_hit == NULL ) return NULL; // SJA:FIXME: should be replaced with an exception
+
+ //gear::Vector3D U(1.0,plane_hit.getU()[1],plane_hit.getU()[0],gear::Vector3D::spherical);
+ //gear::Vector3D V(1.0,plane_hit.getV()[1],plane_hit.getV()[0],gear::Vector3D::spherical);
+ if (type[CEPCConf::TrkHitTypeBit::PLANAR]) {
+ gear::Vector3D U(1.0,trkhit.getCovMatrix(1),trkhit.getCovMatrix(0),gear::Vector3D::spherical);
+ gear::Vector3D V(1.0,trkhit.getCovMatrix(4),trkhit.getCovMatrix(3),gear::Vector3D::spherical);
+ gear::Vector3D X(1.0,0.0,0.0);
+ gear::Vector3D Y(0.0,1.0,0.0);
+ gear::Vector3D Z(0.0,0.0,1.0);
+
+ const float eps = 1.0e-07;
+ // only require vertical to Z axis
+ if (U.dot(Z) > eps) {
+ std::cout << "CEPCArcDiscMeasLayer: TrackerHit measurment vectors U is not vertical to the global Z axis. \n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << std::endl;
+ exit(1);
+ }
+ if (V.dot(Z) > eps) {
+ std::cout << "CEPCArcDiscMeasLayer: TrackerHit measurment vectors V is not vertical to the global Z axis. \n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << std::endl;
+ exit(1);
+ }
+ }
+ /*
+ // U must be the global X axis
+ if( fabs(1.0 - U.dot(Y)) > eps ) {
+ std::cout << "CEPCArcDiscMeasLayer: TrackerHitPlane measurment vectors U is not equal to the global Y axis. \n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << std::endl;
+ exit(1);
+ }
+
+ // V must be the global X axis
+ if( fabs(1.0 - V.dot(X)) > eps ) {
+ std::cout << "CEPCArcDiscMeasLayer: TrackerHitPlane measurment vectors V is not equal to the global X axis. \n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << std::endl;
+ exit(1);
+ }
+ */
+ const edm4hep::Vector3d& pos=trkhit.getPosition();
+ const TVector3 hit(pos.x, pos.y, pos.z);
+
+ // convert to layer coordinates
+ TKalMatrix h = this->XvToMv(hit);
+
+ double x[2] ;
+ double dx[2] ;
+
+ x[0] = h(0, 0);
+ x[1] = h(1, 0);
+
+ if (type[CEPCConf::TrkHitTypeBit::PLANAR]) {
+ dx[0] = trkhit.getCovMatrix(2);
+ dx[1] = trkhit.getCovMatrix(5);
+ }
+ else if (type[CEPCConf::TrkHitTypeBit::CYLINDER]) {
+ dx[0] = sqrt(trkhit.getCovMatrix(0)+trkhit.getCovMatrix(2));
+ dx[1] = sqrt(trkhit.getCovMatrix(5));
+ }
+ else {
+ dx[0] = sqrt(trkhit.getCovMatrix(0)+trkhit.getCovMatrix(2));
+ dx[1] = sqrt(trkhit.getCovMatrix(5));
+ }
+
+ bool hit_on_surface = IsOnSurface(hit);
+#define DEBUG_CONVERT 1
+#ifdef DEBUG_CONVERT
+ std::cout << "CEPCArcDiscMeasLayer::ConvertLCIOTrkHit ILDPlanarHit created"
+ << " u = " << x[0]
+ << " v = " << x[1]
+ << " du = " << dx[0]
+ << " dv = " << dx[1]
+ << " x = " << pos.x
+ << " y = " << pos.y
+ << " z = " << pos.z
+ << " onSurface = " << hit_on_surface
+ << std::endl ;
+#endif
+ return hit_on_surface ? new ILDPlanarHit( *this , x, dx, this->GetBz(), trkhit) : NULL;
+}
diff --git a/Utilities/KalDet/src/ild/common/CEPCCylinderMeasLayer.cc b/Utilities/KalDet/src/ild/common/CEPCCylinderMeasLayer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1dcac2aef163ddf84bae2d1e8b58b98b2caa52ff
--- /dev/null
+++ b/Utilities/KalDet/src/ild/common/CEPCCylinderMeasLayer.cc
@@ -0,0 +1,53 @@
+/** User defined KalTest measurement layer class
+ *
+ * @author
+ */
+
+
+#include "kaltest/TKalTrack.h"
+
+#include "kaldet/CEPCCylinderMeasLayer.h"
+#include "kaldet/ILDCylinderHit.h"
+
+#include <lcio.h>
+#include <edm4hep/TrackerHit.h>
+//#include <EVENT/TrackerHitZCylinder.h>
+
+#include "GaudiKernel/CommonMessaging.h"
+// #include "streamlog/streamlog.h"
+
+#include "TMath.h"
+#include <cmath>
+
+CEPCCylinderMeasLayer::CEPCCylinderMeasLayer(TMaterial &min,
+ TMaterial &mout,
+ Double_t r0,
+ Double_t lhalf,
+ Double_t phi0,
+ Double_t width,
+ Double_t x0,
+ Double_t y0,
+ Double_t z0,
+ Double_t Bz,
+ Bool_t is_active,
+ Int_t CellID,
+ const Char_t *name)
+: ILDCylinderMeasLayer(min, mout, r0, lhalf, x0, y0, z0, Bz, is_active, CellID, name),
+ _phi0(phi0) {
+ _dphi = width/r0/2.0;
+}
+
+Bool_t CEPCCylinderMeasLayer::IsOnSurface(const TVector3 &xx) const {
+
+ bool z = (xx.Z() >= GetZmin() && xx.Z() <= GetZmax());
+ bool r = std::fabs( (xx-this->GetXc()).Perp() - this->GetR() ) < 1.e-3; // for very short, very stiff tracks this can be poorly defined, so we relax this here a bit to 1 micron
+ Double_t phi = (xx-this->GetXc()).Phi();
+ Double_t dphi = phi - _phi0;
+ if (dphi>M_PI) dphi -= M_PI*2;
+ else if (dphi<-M_PI) dphi += M_PI*2;
+ // streamlog_out(DEBUG0) << "CEPCCylinderMeasLayer IsOnSurface for " << this->TVMeasLayer::GetName() << " R = " << this->GetR() << " GetZmin() = " << GetZmin() << " GetZmax() = " << GetZmax()
+ // << " dr = " << std::fabs( (xx-this->GetXc()).Perp() - this->GetR() ) << " r = " << r << " z = " << z
+ // << std::endl;
+
+ return r && z && (fabs(dphi)<_dphi);
+}
diff --git a/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc b/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc
index c9984dfee4cb7845bdd4ef2d7a683c485f79abd5..873803c4cfaf2f4c60f0bc52ca720c0f12dc7e6b 100644
--- a/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc
+++ b/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc
@@ -191,11 +191,13 @@ Bool_t ILDPlanarMeasLayer::IsOnSurface(const TVector3 &xx) const
{
Double_t xi = (xx.Y()-GetXc().Y())*GetNormal().X()/GetNormal().Perp() - (xx.X()-GetXc().X())*GetNormal().Y()/GetNormal().Perp() ;
Double_t zeta = xx.Z() - GetXc().Z();
-
+
bool onSurface = false ;
-
+
if( (xx.X()-GetXc().X())*GetNormal().X() + (xx.Y()-GetXc().Y())*GetNormal().Y() < 1e-4){
- if( xi <= GetXioffset() + GetXiwidth()/2 && xi >= GetXioffset() - GetXiwidth()/2 && TMath::Abs(zeta) <= GetZetawidth()/2){
+ //if( xi <= GetXioffset() + GetXiwidth()/2 && xi >= GetXioffset() - GetXiwidth()/2 && TMath::Abs(zeta) <= GetZetawidth()/2){
+ //FIXME: pick up hits in sensor A but propagate to sensor B
+ if ( xi <= GetXioffset() + GetXiwidth()/2 && xi >= GetXioffset() - GetXiwidth()/2 && TMath::Abs(zeta) <= GetZetawidth()/2 + 10) {
onSurface = true;
}
diff --git a/Utilities/KalDet/src/ild/common/MaterialDataBase.cc b/Utilities/KalDet/src/ild/common/MaterialDataBase.cc
index 50a989e051044936de5d3d51ffaddc43b734ec9e..35281e92cc8ac7ce184e491e22d7b2e9f81b27c6 100644
--- a/Utilities/KalDet/src/ild/common/MaterialDataBase.cc
+++ b/Utilities/KalDet/src/ild/common/MaterialDataBase.cc
@@ -156,7 +156,8 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
TMaterial &beryllium = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.) ;
this->addMaterial(&beryllium, name);
-
+ const double kg_m3TOg_cm3 = 1000.0/1000000.;
+ const double mmTOcm = 1.0/10.;
// TPC Gas
A = 39.948*0.9+(12.011*0.2+1.00794*0.8)*0.1;
Z = 16.4;
@@ -164,6 +165,17 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
// radlen = 1.196e4*2;
radlen = 0.5*1.196e4*2; // SJA:FIXME: reduce by a factor of 10%
name = "tpcgas" ;
+
+ try {
+ const gear::SimpleMaterial& tpc_gas_mat = gearMgr.getSimpleMaterial("TPCGasMaterial");
+ A = tpc_gas_mat.getA();
+ Z = tpc_gas_mat.getZ();
+ density = tpc_gas_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = tpc_gas_mat.getRadLength() * mmTOcm;
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material TPCGasMaterial from GeomSvc! GearMgr=" << &gearMgr << ", will use default gas" << std::endl;
+ }
TMaterial &tpcgas = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&tpcgas, name);
@@ -174,6 +186,17 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
density = 0.0738148 ;
radlen = 489.736 * 0.5 ; // SJA:FIXME just use factor of two for now as the amount differs by ~ factor of 2 from observation in GEANT4
name = "tpcinnerfieldcage" ;
+
+ try {
+ const gear::SimpleMaterial& tpc_inner_mat = gearMgr.getSimpleMaterial("TPCInnerWallMaterial");
+ A = tpc_inner_mat.getA();
+ Z = tpc_inner_mat.getZ();
+ density = tpc_inner_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = tpc_inner_mat.getRadLength() * mmTOcm;
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material TPCInnerWallMaterial from GeomSvc! GearMgr=" << &gearMgr << ", will use default gas" << std::endl;
+ }
TMaterial &tpcinnerfieldcage = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&tpcinnerfieldcage, name);
@@ -184,6 +207,17 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
density = 0.0738148 ;
radlen = 489.736 ;
name = "tpcouterfieldcage" ;
+
+ try {
+ const gear::SimpleMaterial& tpc_outer_mat = gearMgr.getSimpleMaterial("TPCOuterWallMaterial");
+ A = tpc_outer_mat.getA();
+ Z = tpc_outer_mat.getZ();
+ density = tpc_outer_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = tpc_outer_mat.getRadLength() * mmTOcm;
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material TPCOuterWallMaterial from GeomSvc! GearMgr=" << &gearMgr << ", will use default gas" << std::endl;
+ }
TMaterial &tpcouterfieldcage = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&tpcouterfieldcage, name);
@@ -195,8 +229,8 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
A = 12.01 ;
Z = 6.0 ;
- density = 0.5 * 2.00 ; // g/cm^3
- radlen = 2.0 * 21.3485 ; // cm
+ density = /*0.5 */ 2.00 ; // g/cm^3
+ radlen = /*2.0 */ 21.3485 ; // cm
name = "FTDSupportMaterial" ;
TMaterial &ftdsupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.) ;
@@ -216,8 +250,8 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
const gear::SimpleMaterial& vxd_sup_mat = gearMgr.getSimpleMaterial("VXDSupportMaterial");
A = vxd_sup_mat.getA();
Z = vxd_sup_mat.getZ();
- density = vxd_sup_mat.getDensity() * (1000.0/ 1000000.0); // kg/m^3 -> g/cm^3
- radlen = vxd_sup_mat.getRadLength() / 10.0 ; // mm -> cm
+ density = vxd_sup_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = vxd_sup_mat.getRadLength() * mmTOcm;
name = vxd_sup_mat.getName() ;
TMaterial &vxdsupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&vxdsupport, name);
@@ -225,12 +259,25 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
catch( gear::UnknownParameterException& e){
std::cout << "Warning! cannot get material VXDSupportMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
}
+ try {
+ const gear::SimpleMaterial& vxd_bent_mat = gearMgr.getSimpleMaterial("VXDBentSupportMaterial");
+ A = vxd_bent_mat.getA();
+ Z = vxd_bent_mat.getZ();
+ density = vxd_bent_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = vxd_bent_mat.getRadLength() * mmTOcm;
+ name = vxd_bent_mat.getName() ;
+ TMaterial &bentsupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
+ this->addMaterial(&bentsupport, name);
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material OTKBarrelSensorMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
+ }
try {
const gear::SimpleMaterial& vxd_shell_mat = gearMgr.getSimpleMaterial("VXDShellMaterial");
A = vxd_shell_mat.getA();
Z = vxd_shell_mat.getZ();
- density = vxd_shell_mat.getDensity() * (1000.0/ 1000000.0); // kg/m^3 -> g/cm^3
- radlen = vxd_shell_mat.getRadLength() / 10.0 ; // mm -> cm
+ density = vxd_shell_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = vxd_shell_mat.getRadLength() * mmTOcm;
name = vxd_shell_mat.getName() ;
TMaterial &vxdshell = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&vxdshell, name);
@@ -238,12 +285,64 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
catch( gear::UnknownParameterException& e){
std::cout << "Warning! cannot get material VXDShellMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
}
+ try {
+ const gear::SimpleMaterial& itk_sensor_mat = gearMgr.getSimpleMaterial("ITKBarrelSensorMaterial");
+ A = itk_sensor_mat.getA();
+ Z = itk_sensor_mat.getZ();
+ density = itk_sensor_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = itk_sensor_mat.getRadLength() * mmTOcm;
+ name = itk_sensor_mat.getName() ;
+ TMaterial &itksensor = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
+ this->addMaterial(&itksensor, name);
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material ITKBarrelSensorMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
+ }
+ try {
+ const gear::SimpleMaterial& itk_support_mat = gearMgr.getSimpleMaterial("ITKBarrelSupportMaterial");
+ A = itk_support_mat.getA();
+ Z = itk_support_mat.getZ();
+ density = itk_support_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = itk_support_mat.getRadLength() * mmTOcm;
+ name = itk_support_mat.getName() ;
+ TMaterial &itksupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
+ this->addMaterial(&itksupport, name);
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material ITKBarrelSupportMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
+ }
+ try {
+ const gear::SimpleMaterial& itk_support_mat = gearMgr.getSimpleMaterial("ITKEndcapSupportMaterial");
+ A = itk_support_mat.getA();
+ Z = itk_support_mat.getZ();
+ density = itk_support_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = itk_support_mat.getRadLength() * mmTOcm;
+ name = itk_support_mat.getName() ;
+ TMaterial &itksupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
+ this->addMaterial(&itksupport, name);
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material ITKEndcapSupportMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
+ }
+ try {
+ const gear::SimpleMaterial& otk_support_mat = gearMgr.getSimpleMaterial("OTKBarrelSupportMaterial");
+ A = otk_support_mat.getA();
+ Z = otk_support_mat.getZ();
+ density = otk_support_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = otk_support_mat.getRadLength() * mmTOcm;
+ name = otk_support_mat.getName() ;
+ TMaterial &otksupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
+ this->addMaterial(&otksupport, name);
+ }
+ catch( gear::UnknownParameterException& e){
+ std::cout << "Warning! cannot get material OTKBarrelSensorMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl;
+ }
try {
const gear::SimpleMaterial& etd_sup_mat = gearMgr.getSimpleMaterial("OTKEndcapSupportMaterial");
A = etd_sup_mat.getA();
Z = etd_sup_mat.getZ();
- density = etd_sup_mat.getDensity() * (1000.0/ 1000000.0); // kg/m^3 -> g/cm^3
- radlen = etd_sup_mat.getRadLength() / 10.0 ; // mm -> cm
+ density = etd_sup_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = etd_sup_mat.getRadLength() * mmTOcm;
name = etd_sup_mat.getName() ;
TMaterial &etdsupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&etdsupport, name);
@@ -255,8 +354,8 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
const gear::SimpleMaterial& tpc_readout_mat = gearMgr.getSimpleMaterial("TPCReadoutMaterial");
A = tpc_readout_mat.getA();
Z = tpc_readout_mat.getZ();
- density = tpc_readout_mat.getDensity() * (1000.0/ 1000000.0); // kg/m^3 -> g/cm^3
- radlen = tpc_readout_mat.getRadLength() / 10.0 ; // mm -> cm
+ density = tpc_readout_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = tpc_readout_mat.getRadLength() * mmTOcm;
name = tpc_readout_mat.getName() ;
TMaterial &tpcreadout = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&tpcreadout, name);
@@ -268,8 +367,8 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
const gear::SimpleMaterial& tpc_endplate_mat = gearMgr.getSimpleMaterial("TPCEndplateMaterial");
A = tpc_endplate_mat.getA();
Z = tpc_endplate_mat.getZ();
- density = tpc_endplate_mat.getDensity() * (1000.0/ 1000000.0); // kg/m^3 -> g/cm^3
- radlen = tpc_endplate_mat.getRadLength() / 10.0 ; // mm -> cm
+ density = tpc_endplate_mat.getDensity() * kg_m3TOg_cm3;
+ radlen = tpc_endplate_mat.getRadLength() * mmTOcm;
name = tpc_endplate_mat.getName() ;
TMaterial &tpcendplate = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.);
this->addMaterial(&tpcendplate, name);
diff --git a/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc b/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc
index ffe84e4433f2f9608e04a11f9825b978987b1069..b2848bd24b2405f41726cc2773276e63f5536b5a 100644
--- a/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc
+++ b/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc
@@ -39,11 +39,11 @@ CEPCOTKKalDetector::CEPCOTKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
TMaterial & air = *MaterialDataBase::Instance().getMaterial("air");
TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon");
- TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon");
+ TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("OTKBarrelSupportMaterial");
if(geoSvc){
setupGearGeom( geoSvc ) ;
- }
+ }
else{
setupGearGeom( gearMgr );
}
diff --git a/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc b/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc
index e4a6b179dc85c3a2b9ba9a681f4f087c15a5785e..e85f3ba1fa3d478beabcb874d0f0b80415f5d478 100644
--- a/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc
+++ b/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc
@@ -33,8 +33,10 @@ CEPCITKKalDetector::CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
MaterialDataBase::Instance().registerForService(gearMgr, geoSvc);
TMaterial & air = *MaterialDataBase::Instance().getMaterial("air");
- TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon");
- TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon");
+ //TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon");
+ TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("ITKBarrelSensorMaterial");
+ //TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon");
+ TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("ITKBarrelSupportMaterial");
if(geoSvc){
this->setupGearGeom(geoSvc);
@@ -70,13 +72,13 @@ CEPCITKKalDetector::CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
const double phi0 = _ITKgeo[layer].phi0 ;
- const double ladder_distance = _ITKgeo[layer].supRMin ;
+ const double ladder_distance = _ITKgeo[layer].supRMin ;
const double ladder_thickness = _ITKgeo[layer].supThickness ;
- const double sensitive_distance = _ITKgeo[layer].senRMin ;
+ const double sensitive_distance = _ITKgeo[layer].senRMin ;
const double sensitive_thickness = _ITKgeo[layer].senThickness ;
- const double width = _ITKgeo[layer].width ;
+ const double width = _ITKgeo[layer].width ;
const double length = _ITKgeo[layer].length;
double currPhi;
@@ -87,7 +89,9 @@ CEPCITKKalDetector::CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
const int nsensors = _ITKgeo[layer].nSensorsPerLadder;
const double sensor_length = _ITKgeo[layer].sensorLength;
-
+
+ const int nrow = std::floor((sensor_length*nsensors)/length+0.001);
+ const double gap = length - sensor_length*nsensors/nrow;
// TODO: sorting_policy for overlap region like VXD
for (int ladder=0; ladder<nLadders; ++ladder) {
@@ -117,7 +121,8 @@ CEPCITKKalDetector::CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+1) * eps_layer + eps_sensor * isensor ;
- double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor*sensor_length) ;
+ double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor%(nsensors/nrow))*sensor_length;
+ if (z_centre_sensor>0) z_centre_sensor += gap;
if (_isStripDetector) {
// measurement plane defined as the middle of the sensitive volume
@@ -156,7 +161,9 @@ CEPCITKKalDetector::CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+2) * eps_layer + eps_sensor * isensor ;
- double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor*sensor_length) ;
+ //double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor*sensor_length) ;
+ double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor%(nsensors/nrow))*sensor_length;
+ if (z_centre_sensor>0) z_centre_sensor += gap;
if (_isStripDetector) {
// measurement plane defined as the middle of the sensitive volume
@@ -242,9 +249,8 @@ void CEPCITKKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){
_ITKgeo[layer].nSensorsPerLadder = 1 ;
}
- _ITKgeo[layer].sensorLength = _ITKgeo[layer].length / _ITKgeo[layer].nSensorsPerLadder;
-
-
+ _ITKgeo[layer].sensorLength = pITKDetMain.getDoubleVals("length_sensors")[layer];//_ITKgeo[layer].length / _ITKgeo[layer].nSensorsPerLadder;
+
if (_isStripDetector) {
_ITKgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ;
} else {
@@ -303,7 +309,7 @@ void CEPCITKKalDetector::setupGearGeom( IGeomSvc* geoSvc ){
_ITKgeo[layer].senThickness = pITKLayerLayout.thicknessSensitive*CLHEP::cm;
_ITKgeo[layer].supThickness = pITKLayerLayout.thicknessSupport*CLHEP::cm;
_ITKgeo[layer].nSensorsPerLadder = pITKLayerLayout.sensorsPerLadder;
- _ITKgeo[layer].sensorLength = _ITKgeo[layer].length / _ITKgeo[layer].nSensorsPerLadder;
+ _ITKgeo[layer].sensorLength = pITKLayerLayout.lengthSensor*CLHEP::cm;//_ITKgeo[layer].length / _ITKgeo[layer].nSensorsPerLadder;
if (_isStripDetector) {
_ITKgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ;
diff --git a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc
index 44484c7dc7bb550712f55f0032b10ccf69e6ff9b..48b73b4148102fd3ed2c2e3e4cd63e608502f796 100644
--- a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc
+++ b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc
@@ -330,8 +330,8 @@ TVKalDetector(10)
double r_max_ecal_ecap = ecalE.getExtent()[1];
double z_min_ecal_ecap = ecalE.getExtent()[2];
+ encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::ECAL_ENDCAP ;
encoder[lcio::ILDCellID0::module] = 0;
-
encoder[lcio::ILDCellID0::side] = lcio::ILDDetID::fwd;
TVector3 front_face_centre_fwd( 0.0, 0.0, z_min_ecal_ecap); // for +z
diff --git a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc
index 8f9a33043186d5878b583954942f35c1c1a7b33e..847838ecd9a6c7187c0c7126202c211c5c2c9d40 100644
--- a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc
+++ b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc
@@ -3,7 +3,7 @@
#include "kaldet/MaterialDataBase.h"
#include "kaldet/ILDParallelPlanarMeasLayer.h"
-#include "kaldet/ILDCylinderMeasLayer.h"
+#include "kaldet/CEPCCylinderMeasLayer.h"
#include "kaldet/ILDDiscMeasLayer.h"
#include <UTIL/BitField64.h>
@@ -211,9 +211,12 @@ CEPCVTXKalDetector::CEPCVTXKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
}
}
- for (int layer=0; layer<_nLayers[1]; ++layer) {
+ if (_nLayers[1]>0) {
+ TMaterial & bentmat = *MaterialDataBase::Instance().getMaterial("VXDBentSupportMaterial");
+ for (int layer=0; layer<_nLayers[1]; ++layer) {
//std::cout << "add stitching ... " << layer << std::endl;
double phi0 = _STTgeo[layer].phi0;
+ double width = _STTgeo[layer].width;
double ladder_distance = _STTgeo[layer].supRMin;
double ladder_thickness = _STTgeo[layer].supThickness;
@@ -249,16 +252,17 @@ CEPCVTXKalDetector::CEPCVTXKalDetector( const gear::GearMgr& gearMgr, IGeomSvc*
}
// air - sensitive boundary
- Add(new ILDCylinderMeasLayer(air, silicon, sensitive_distance, halfz, x0, y0, z0, _bZ, dummy, -1, "STTMeasL_0"));
+ Add(new CEPCCylinderMeasLayer(air, silicon, sensitive_distance, halfz, currPhi, width, 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"));
+ Add(new CEPCCylinderMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*(1.-_relative_position_of_measurement_surface),
+ halfz, currPhi, width, 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"));
+ Add(new CEPCCylinderMeasLayer(silicon, bentmat, sensitive_distance+sensitive_thickness, halfz, currPhi, width, 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" )) ;
+ Add(new CEPCCylinderMeasLayer(bentmat, air, ladder_distance+ladder_thickness, halfz, currPhi, width, x0, y0, z0, _bZ, dummy, -1, "STTSupRear_0" )) ;
}
+ }
}
SetOwner();
@@ -307,12 +311,14 @@ void CEPCVTXKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){
const std::vector<double> phi0s = pVXDDetMain.getDoubleVals("VTXLayerPhi0");
const std::vector<double> rgaps = pVXDDetMain.getDoubleVals("VTXLayerRadialGap");
const std::vector<double> dphis = pVXDDetMain.getDoubleVals("VTXLayerDeltaPhi");
+ const std::vector<double> widths = pVXDDetMain.getDoubleVals("VTXLayerWidth");
_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].width = widths[ilayer];
_STTgeo[ilayer].senRMin = rsens[ilayer];
_STTgeo[ilayer].senThickness = tsens[ilayer];
_STTgeo[ilayer].supRMin = rsups[ilayer];