diff --git a/Detector/DetCRD/CMakeLists.txt b/Detector/DetCRD/CMakeLists.txt index d3f0b06954b5d3a547cb13495dfb52201ebf0818..de22c5023717faaafc4938d5bc6876b780cf8261 100644 --- a/Detector/DetCRD/CMakeLists.txt +++ b/Detector/DetCRD/CMakeLists.txt @@ -21,6 +21,7 @@ gaudi_add_module(DetCRD src/Calorimeter/RotatedCrystalCalorimeter_v01_geo.cpp src/Other/CRDBeamPipe_v01_geo.cpp src/Tracker/SiTrackerSkewRing_v01_geo.cpp + src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp LINK ${DD4hep_COMPONENT_LIBRARIES} ) diff --git a/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_01.xml index 9f9aa7edc209ad3be8bcffea20edb0399dac3303..1ef2c4b96676946979cddee80bb62df06356e012 100644 --- a/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_01.xml +++ b/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_01.xml @@ -13,53 +13,54 @@ </define> <detectors> - <detector name="BeamPipe" type="DD4hep_CRDBeamPipe_v01" vis="BeamPipeVis"> + <detector name="BeamPipe" type="CRDBeamPipe_v01" vis="VacVis"> <parameter crossingangle="CrossingAngle" /> <envelope> <shape type="Assembly"/> </envelope> <section type ="Center" name="IPInnerTube" zStart="0" zEnd="BeamPipe_CentralBe_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Central_inner_radius"/> - <layer material="G4_Be" thickness="BeamPipe_Be_inner_thickness"/> - <layer material="G4_PARAFFIN" thickness="BeamPipe_Cooling_thickness"/> - <layer material="G4_Be" thickness="BeamPipe_Be_outer_thickness"/> + <layer material="beam" thickness="BeamPipe_Central_inner_radius" vis="VacVis"/> + <layer material="G4_Be" thickness="BeamPipe_Be_inner_thickness" vis="TubeVis"/> + <layer material="G4_PARAFFIN" thickness="BeamPipe_Cooling_thickness" vis="GrayVis"/> + <layer material="G4_Be" thickness="BeamPipe_Be_outer_thickness" vis="TubeVis"/> </section> <section type="Center" name="IPAl" zStart="BeamPipe_CentralBe_zmax" zEnd="BeamPipe_CentralAl_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Central_inner_radius"/> - <layer material="G4_Al" thickness="BeamPipe_Al_thickness"/> + <layer material="beam" thickness="BeamPipe_Central_inner_radius" vis="VacVis"/> + <layer material="G4_Al" thickness="BeamPipe_Al_thickness" vis="TubeVis"/> </section> <section type="Center" name="ExpandPipe" zStart="BeamPipe_CentralAl_zmax" zEnd="BeamPipe_ConeAl_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Central_inner_radius" thicknessEnd="BeamPipe_Expanded_inner_radius"/> - <layer material="G4_Al" thickness="BeamPipe_Al_thickness" thicknessEnd="BeamPipe_Al_thickness"/> + <layer material="beam" thickness="BeamPipe_Central_inner_radius" thicknessEnd="BeamPipe_Expanded_inner_radius" vis="VacVis"/> + <layer material="G4_Al" thickness="BeamPipe_Al_thickness" thicknessEnd="BeamPipe_Al_thickness" vis="TubeVis"/> </section> <section type="Center" name="ThickPipe" zStart="BeamPipe_ConeAl_zmax" zEnd="BeamPipe_LinkerAl_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Expanded_inner_radius"/> - <layer material="G4_Al" thickness="BeamPipe_Al_thickness"/> + <layer material="beam" thickness="BeamPipe_Expanded_inner_radius" vis="VacVis"/> + <layer material="G4_Al" thickness="BeamPipe_Al_thickness" vis="TubeVis"/> </section> <section type="CenterSide" name="OutsideLink" zStart="BeamPipe_LinkerAl_zmax" zEnd="BeamPipe_LinkerCu_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Expanded_inner_radius"/> - <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness"/> + <layer material="beam" thickness="BeamPipe_Expanded_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> </section> <section type="FatWaist" name="Waist" zStart="BeamPipe_LinkerCu_zmax" zEnd="BeamPipe_Waist_zmax" rStart="BeamPipe_Expanded_inner_radius" size="BeamPipe_Crotch_hole_height"> - <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> </section> - <!--CrotchAsymUp&CrotchAsymDn not work to fix, because of problem on convert from TGeo to Geant4--> - <!--section type="CrotchAsymUp" name="Fork" zStart="BeamPipe_Waist_zmax" zEnd="BeamPipe_Crotch_zmax" + <!--CrotchAsymUp&CrotchAsymDn not work to fix, because of problem on convert from TGeo to Geant4--> + <!--Since lcg101, they work--> + <section type="CrotchAsymUp" name="Fork" zStart="BeamPipe_Waist_zmax" zEnd="BeamPipe_Crotch_zmax" rStart="BeamPipe_Expanded_inner_radius" rEnd="BeamPipe_Upstream_inner_radius" size="BeamPipe_Crotch_hole_height"> - <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" thicknessEnd="ForkAsymThickness"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" thicknessEnd="ForkAsymThickness" vis="TubeVis"/> </section> <section type="CrotchAsymDn" name="Fork" zStart="BeamPipe_Waist_zmax" zEnd="BeamPipe_Crotch_zmax" rStart="BeamPipe_Expanded_inner_radius" rEnd="BeamPipe_Dnstream_inner_radius" size="BeamPipe_Crotch_hole_height"> - <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness"/> - </section--> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> <section type="FlareLegUp" name="FirstDoublePipe" zStart="BeamPipe_Crotch_zmax" zEnd="BeamPipe_FirstSeparated_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Upstream_inner_radius" thicknessEnd="BeamPipe_Dnstream_inner_radius"/> - <layer material="G4_Cu" thickness="ForkAsymThickness" thicknessEnd="BeamPipe_Cu_thickness"/> + <layer material="beam" thickness="BeamPipe_Upstream_inner_radius" thicknessEnd="BeamPipe_Dnstream_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="ForkAsymThickness" thicknessEnd="BeamPipe_Cu_thickness" vis="TubeVis"/> </section> <section type="FlareLegDn" name="FirstDoublePipe" zStart="BeamPipe_Crotch_zmax" zEnd="BeamPipe_FirstSeparated_zmax" rStart="0"> - <layer material="beam" thickness="BeamPipe_Dnstream_inner_radius"/> - <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness"/> + <layer material="beam" thickness="BeamPipe_Dnstream_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> </section> </detector> </detectors> diff --git a/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_02.xml b/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_02.xml new file mode 100644 index 0000000000000000000000000000000000000000..ef58b72abe1e9d6e836043f4b03a7cfb0d2316bc --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v01/Beampipe_v01_02.xml @@ -0,0 +1,68 @@ +<lccdd> + <info name="CRD" title="CRD Beam pipe" author="Chengdong Fu" url="no" status="development" version="1.0"> + <comment>A beampipe for CRD</comment> + </info> + + <display> + <vis name="TubeVis" alpha="0.1" r="1.0" g="0.7" b="0.5" showDaughters="true" visible="true"/> + <vis name="VacVis" alpha="1.0" r="0.0" g="0.0" b="0.0" showDaughters="true" visible="false"/> + </display> + + <define> + <!--only needed for asymetry double pipe--> + <!--constant name="ForkAsymThickness" value="BeamPipe_Dnstream_inner_radius+BeamPipe_Cu_thickness-BeamPipe_Upstream_inner_radius"/--> + </define> + + <detectors> + <detector name="BeamPipe" type="CRDBeamPipe_v01" vis="VacVis"> + <parameter crossingangle="CrossingAngle" /> + <envelope> + <shape type="Assembly"/> + </envelope> + + <section type ="Center" name="IPInnerTube" zStart="0" zEnd="BeamPipe_CentralBe_zmax" rStart="0"> + <layer material="beam" thickness="BeamPipe_Central_inner_radius" vis="VacVis"/> + <layer material="G4_Be" thickness="BeamPipe_Be_inner_thickness" vis="TubeVis"/> + <layer material="G4_PARAFFIN" thickness="BeamPipe_Cooling_thickness" vis="GrayVis"/> + <layer material="G4_Be" thickness="BeamPipe_Be_outer_thickness" vis="TubeVis"/> + </section> + <section type="Center" name="IPAl" zStart="BeamPipe_CentralBe_zmax" zEnd="BeamPipe_CentralAl_zmax" rStart="0"> + <layer material="beam" thickness="BeamPipe_Central_inner_radius" vis="VacVis"/> + <layer material="G4_Al" thickness="BeamPipe_Al_thickness" vis="TubeVis"/> + </section> + <section type="Waist" name="Waist1st" zStart="BeamPipe_CentralAl_zmax" zEnd="BeamPipe_ExpandAl_zmax" rStart="BeamPipe_Central_inner_radius" size="BeamPipe_FirstExpand_width"> + <layer material="G4_Al" thickness="BeamPipe_Al_thickness" vis="TubeVis"/> + </section> + <section type="Runway" name="Waist2nd" zStart="BeamPipe_ExpandAl_zmax" zEnd="BeamPipe_Linker_zmin" rStart="BeamPipe_Central_inner_radius" size="BeamPipe_FirstExpand_width"> + <layer material="G4_Al" thickness="BeamPipe_Al_thickness" vis="TubeVis"/> + </section> + <section type="Runway" name="Waist3rd" zStart="BeamPipe_Linker_zmin" zEnd="BeamPipe_Linker_zmax" rStart="BeamPipe_Central_inner_radius" size="BeamPipe_FirstExpand_width"> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> + <section type="Runway" name="Waist4th" zStart="BeamPipe_Linker_zmax" zEnd="BeamPipe_Waist_zmax" rStart="BeamPipe_Central_inner_radius" size="BeamPipe_FirstExpand_width" + shift="BeamPipe_SecondExpand_width-BeamPipe_FirstExpand_width"> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> + <section type="Crotch" name="Fork" zStart="BeamPipe_Waist_zmax" zEnd="BeamPipe_Crotch_zmax" + rStart="BeamPipe_Central_inner_radius" rEnd="BeamPipe_Central_inner_radius" size="BeamPipe_SecondExpand_width"> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> + <section type="Legs" name="FirstDoublePipe" zStart="BeamPipe_Crotch_zmax" zEnd="BeamPipe_FirstSeparated_zmax" rStart="0"> + <layer material="beam" thickness="BeamPipe_Fork_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> + <section type="Legs" name="BeforeMask" zStart="BeamPipe_FirstSeparated_zmax" zEnd="BeamPipe_Mask_zmin" rStart="0"> + <layer material="beam" thickness="BeamPipe_Fork_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> + <section type="Legs" name="Mask" zStart="BeamPipe_Mask_zmin" zEnd="BeamPipe_Mask_zmax" rStart="0"> + <layer material="beam" thickness="BeamPipe_Mask_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness+BeamPipe_Fork_inner_radius-BeamPipe_Mask_inner_radius" vis="TubeVis"/> + </section> + <section type="Legs" name="SecondDoublePipe" zStart="BeamPipe_Mask_zmax" zEnd="BeamPipe_SecondSeparated_zmax" rStart="0"> + <layer material="beam" thickness="BeamPipe_Fork_inner_radius" vis="VacVis"/> + <layer material="G4_Cu" thickness="BeamPipe_Cu_thickness" vis="TubeVis"/> + </section> + </detector> + </detectors> +</lccdd> diff --git a/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml new file mode 100644 index 0000000000000000000000000000000000000000..2c3c63d33dfebfe092b71deddd759918a9a39228 --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml @@ -0,0 +1,107 @@ +<lccdd> + <info name="VXD_StaggeredLadder_v01_01" + title="CepC VXD with staggered ladders" + author="H.Zeng, " + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC vertex detector based on MOST2 project </comment> + </info> + <define> + <constant name="VXD_inner_radius" value="Vertex_inner_radius"/> + <constant name="VXD_outer_radius" value="Vertex_outer_radius"/> + <constant name="VXD_half_length" value="Vertex_half_length"/> + <constant name="VXDLayer1_half_length" value="90*mm" /> + <constant name="VXDLayer2_half_length" value="90*mm" /> + <constant name="VXDLayer3_half_length" value="90*mm" /> + <constant name="VXD_sensor_length" value="30*mm" /> + </define> + + <detectors> + <detector id="DetID_VXD" name="VXD" type="SiTrackerStaggeredLadder_v01" vis="VXDVis" readout="VXDCollection" insideTrackingVolume="true"> + <envelope> + <shape type="BooleanShape" operation="Subtraction" material="Air" > + <shape type="BooleanShape" operation="Subtraction" material="Air" > + <shape type="Tube" rmin="VXD_inner_radius+1*mm" rmax="VXD_outer_radius" dz="VXD_half_length" /> + <shape type="Cone" rmin1="0" rmax1="BeamPipe_VertexRegion_rmax" rmin2="0" rmax2="Vertex_Side_rmin" z="(VXD_half_length-BeamPipe_CentralAl_zmax)/2." /> + <position x="0" y="0" z="VXD_half_length-(VXD_half_length-BeamPipe_CentralAl_zmax)/2."/> + </shape> + <shape type="Cone" rmin1="0" rmax1="BeamPipe_VertexRegion_rmax" rmin2="0" rmax2="Vertex_Side_rmin" z="(VXD_half_length-BeamPipe_CentralAl_zmax)/2." /> + <position x="0" y="0" z="-(VXD_half_length-(VXD_half_length-BeamPipe_CentralAl_zmax)/2.)"/> + <rotation x="0" y="180.*deg" z="0" /> + </shape> + </envelope> + + <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_PIXEL "/> + + <global sensitive_thickness="VXD_sensitive_thickness" support_thickness="VXD_support_thickness" sensor_length="VXD_sensor_length" + sensitive_mat="G4_Si" support_mat="G4_C" sensitive_threshold_KeV="64*keV" /> + <display ladder="SeeThrough" support="VXDSupportVis" flex="VXDFlexVis" sens_env="SeeThrough" sens="GrayVis" deadsensor="GreenVis" deadwire="RedVis"/> + + <layer layer_id="0" ladder_radius="17.4*mm" ladder_offset="(8.4-1.5)*mm" n_sensors_per_side="VXDLayer1_half_length*2/VXD_sensor_length" + n_ladders="10" > + <ladder isDoubleSided="true"> + <ladderSupport height="2*mm" length="200*mm" thickness="350*um" width="16.8*mm" mat="CarbonFiber"/> + <flex n_slices="3"> + <slice length="200*mm" thickness="60*um" width="16.8*mm" mat="Epoxy"/> + <slice length="200*mm" thickness="74*um" width="16.8*mm" mat="Kapton"/> + <slice length="200*mm" thickness="26.8*um" width="16.8*mm" mat="G4_Al"/> + </flex> + <sensor n_sensors="7" gap="0.1*mm" thickness="50*um" active_length="25.6*mm" active_width="12.8*mm" dead_width="2*mm" sensor_mat="G4_Si" + deadwire_length="(7*(25.6+0.1)-0.1)*mm" deadwire_width="2*mm" deadwire_thickness="(50/10)*um" deadwire_mat="G4_Al"/> + </ladder> + </layer> + <layer layer_id="1" ladder_radius="36.9*mm" ladder_offset="(8.4+5.0)*mm" n_sensors_per_side="VXDLayer2_half_length*2/VXD_sensor_length" + n_ladders="22" > + <ladder isDoubleSided="true"> + <ladderSupport height="2*mm" length="200*mm" thickness="350*um" width="16.8*mm" mat="CarbonFiber"/> + <flex n_slices="3"> + <slice length="200*mm" thickness="60*um" width="16.8*mm" mat="Epoxy"/> + <slice length="200*mm" thickness="74*um" width="16.8*mm" mat="Kapton"/> + <slice length="200*mm" thickness="26.8*um" width="16.8*mm" mat="G4_Al"/> + <!-- <slice length="200*mm" thickness="15*um" width="16.8*mm" mat="Epoxy"/> --> + <!-- <slice length="200*mm" thickness="12*um" width="16.8*mm" mat="Kapton"/> --> + <!-- <slice length="200*mm" thickness="15*um" width="16.8*mm" mat="Epoxy"/> --> + <!-- <slice length="200*mm" thickness="13.4*um" width="16.8*mm" mat="G4_Al"/> --> + <!-- <slice length="200*mm" thickness="50*um" width="16.8*mm" mat="Kapton"/> --> + <!-- <slice length="200*mm" thickness="13.4*um" width="16.8*mm" mat="G4_Al"/> --> + <!-- <slice length="200*mm" thickness="15*um" width="16.8*mm" mat="Epoxy"/> --> + <!-- <slice length="200*mm" thickness="12*um" width="16.8*mm" mat="Kapton"/> --> + <!-- <slice length="200*mm" thickness="15*um" width="14.8*mm" mat="Epoxy"/> --> + </flex> + <sensor n_sensors="7" gap="0.1*mm" thickness="50*um" active_length="25.6*mm" active_width="12.8*mm" dead_width="2*mm" sensor_mat="G4_Si" + deadwire_length="(7*(25.6+0.1)-0.1)*mm" deadwire_width="2*mm" deadwire_thickness="(50/10)*um" deadwire_mat="G4_Al"/> + </ladder> + </layer> + <layer layer_id="2" ladder_radius="57.7*mm" ladder_offset="(8.4+9.6)*mm" n_sensors_per_side="VXDLayer3_half_length*2/VXD_sensor_length" + n_ladders="32" > + <ladder isDoubleSided="true"> + <ladderSupport height="2*mm" length="200*mm" thickness="350*um" width="16.8*mm" mat="CarbonFiber"/> + <flex n_slices="3"> + <slice length="200*mm" thickness="60*um" width="16.8*mm" mat="Epoxy"/> + <slice length="200*mm" thickness="74*um" width="16.8*mm" mat="Kapton"/> + <slice length="200*mm" thickness="26.8*um" width="16.8*mm" mat="G4_Al"/> + <!-- <slice length="200*mm" thickness="15*um" width="16.8*mm" mat="Epoxy"/> + <slice length="200*mm" thickness="12*um" width="16.8*mm" mat="Kapton"/> + <slice length="200*mm" thickness="15*um" width="16.8*mm" mat="Epoxy"/> + <slice length="200*mm" thickness="13.4*um" width="16.8*mm" mat="G4_Al"/> + <slice length="200*mm" thickness="50*um" width="16.8*mm" mat="Kapton"/> + <slice length="200*mm" thickness="13.4*um" width="16.8*mm" mat="G4_Al"/> + <slice length="200*mm" thickness="15*um" width="16.8*mm" mat="Epoxy"/> + <slice length="200*mm" thickness="12*um" width="16.8*mm" mat="Kapton"/> + <slice length="200*mm" thickness="15*um" width="14.8*mm" mat="Epoxy"/> --> + </flex> + <sensor n_sensors="7" gap="0.1*mm" thickness="50*um" active_length="25.6*mm" active_width="12.8*mm" dead_width="2*mm" sensor_mat="G4_Si" + deadwire_length="(7*(25.6+0.1)-0.1)*mm" deadwire_width="2*mm" deadwire_thickness="(50/10)*um" deadwire_mat="G4_Al"/> + </ladder> + </layer> + </detector> + + </detectors> + + <readouts> + <readout name="VXDCollection"> + <id>system:5,side:-2,layer:9,module:8,active:8,sensor:8</id> + </readout> + </readouts> +</lccdd> diff --git a/Detector/DetCRD/compact/CRD_o1_v03/CRD_o1_v03-onlyVXD.xml b/Detector/DetCRD/compact/CRD_o1_v03/CRD_o1_v03-onlyVXD.xml new file mode 100644 index 0000000000000000000000000000000000000000..44eb691c033035b07cf8d0a435fe8defdeac4014 --- /dev/null +++ b/Detector/DetCRD/compact/CRD_o1_v03/CRD_o1_v03-onlyVXD.xml @@ -0,0 +1,50 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + <info name="CRD_o1_v02" + title="CepC reference detctor with coil inside Hcal, pixel SIT and strip SET" + author="Hao Zeng" + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC reference detector simulation models used for detector study </comment> + </info> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="../CRD_common_v01/materials.xml"/> + <gdmlFile ref="../materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="25*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> + + <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/> + </define> + + <include ref="../CRD_o1_v02/CRD_Dimensions_v01_02.xml"/> + + <include ref="../CRD_common_v01/VXD_StaggeredLadder_v01_01.xml"/> + + <fields> + <field name="InnerSolenoid" type="solenoid" + inner_field="Field_nominal_value" + outer_field="0" + zmax="SolenoidCoil_half_length" + inner_radius="SolenoidCoil_center_radius" + outer_radius="Solenoid_outer_radius"> + </field> + <field name="OuterSolenoid" type="solenoid" + inner_field="0" + outer_field="Field_outer_nominal_value" + zmax="SolenoidCoil_half_length" + inner_radius="Solenoid_outer_radius" + outer_radius="Yoke_barrel_inner_radius"> + </field> + </fields> + +</lccdd> diff --git a/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py b/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py new file mode 100644 index 0000000000000000000000000000000000000000..5d595621af36d9284d9ac4e1f78adf85c0838f5d --- /dev/null +++ b/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python +from Gaudi.Configuration import * + +from Configurables import k4DataSvc +dsvc = k4DataSvc("EventDataSvc") + +from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ + +seed = [10] +# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi +rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4 +rndmengine.SetSingleton = True +rndmengine.Seeds = seed + +rndmgensvc = RndmGenSvc("RndmGenSvc") +rndmgensvc.Engine = rndmengine.name() + +#geometry_option = "CRD_o1_v01/CRD_o1_v01.xml" +geometry_option = "CRD_o1_v03/CRD_o1_v03-onlyVXD.xml" +#... + +if not os.getenv("DETCRDROOT"): + print("Can't find the geometry. Please setup envvar DETCRDROOT." ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option) +if not os.path.exists(geometry_path): + print("Can't find the compact geometry file: %s"%geometry_path) + sys.exit(-1) + +from Configurables import GeomSvc +geosvc = GeomSvc("GeomSvc") +geosvc.compact = geometry_path + +############################################################################## +# Physics Generator +############################################################################## +from Configurables import GenAlgo +from Configurables import GtGunTool +from Configurables import StdHepRdr +from Configurables import SLCIORdr +from Configurables import HepMCRdr +from Configurables import GenPrinter +gun = GtGunTool("GtGunTool") +gun.Particles = ["mu-"] +#gun.Particles = ["nu_e"] +#gun.PositionXs = [0] +#gun.PositionYs = [0] +#gun.PositionZs = [0] +gun.EnergyMins = [100.] # GeV +gun.EnergyMaxs = [100.] # GeV +gun.ThetaMins = [0] # deg +gun.ThetaMaxs = [180] # deg +gun.PhiMins = [0] # deg +gun.PhiMaxs = [360] # deg +# stdheprdr = StdHepRdr("StdHepRdr") +# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" +# lciordr = SLCIORdr("SLCIORdr") +# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio" +# hepmcrdr = HepMCRdr("HepMCRdr") +# hepmcrdr.Input = "example_UsingIterators.txt" + +genprinter = GenPrinter("GenPrinter") + +genalg = GenAlgo("GenAlgo") +genalg.GenTools = ["GtGunTool"] +#genalg.GenTools = ["StdHepRdr"] +# genalg.GenTools = ["StdHepRdr", "GenPrinter"] +# genalg.GenTools = ["SLCIORdr", "GenPrinter"] +# genalg.GenTools = ["HepMCRdr", "GenPrinter"] + +############################################################################## +# Detector Simulation +############################################################################## +from Configurables import DetSimSvc +detsimsvc = DetSimSvc("DetSimSvc") + +from Configurables import DetSimAlg +detsimalg = DetSimAlg("DetSimAlg") +detsimalg.RandomSeeds = seed +# detsimalg.VisMacs = ["vis.mac"] +detsimalg.RunCmds = [ +# "/tracking/verbose 1", +] +detsimalg.AnaElems = [ + # example_anatool.name() +# "ExampleAnaElemTool", + "Edm4hepWriterAnaElemTool" +] +detsimalg.RootDetElem = "WorldDetElemTool" + +# output +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "CRD-skewVXD-Sim-100.root" +out.outputCommands = ["keep *"] + +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [genalg, detsimalg, out], + EvtSel = 'NONE', + EvtMax = 500, + ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc], + OutputLevel=INFO +) diff --git a/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py b/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py index a424a1b3af88ed4b9929cf87ade65891277ce763..1386affd7623cda211e7d1d8f983fa749b26b9f3 100644 --- a/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py +++ b/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py @@ -88,6 +88,10 @@ detsimalg.AnaElems = [ ] detsimalg.RootDetElem = "WorldDetElemTool" +from Configurables import CalorimeterSensDetTool +cal_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool") +cal_sensdettool.CalNamesApplyBirks = ["EcalBarrel"] + # output from Configurables import PodioOutput out = PodioOutput("outputalg") diff --git a/Detector/DetCRD/src/Calorimeter/CRDEcal_v01.cpp b/Detector/DetCRD/src/Calorimeter/CRDEcal_v01.cpp index f75648822cdf44352cd3cec502d292c1157b7d2c..4ba45b249f46c9bd815ac60b21f3a32f385a51cf 100644 --- a/Detector/DetCRD/src/Calorimeter/CRDEcal_v01.cpp +++ b/Detector/DetCRD/src/Calorimeter/CRDEcal_v01.cpp @@ -182,7 +182,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, DetElement sd(ECAL, _toString(i,"trap%3d"), detid); sd.setPlacement(plv); } - + sens.setType("calorimeter"); ECAL.addExtension< LayeredCalorimeterData >( caloData ) ; diff --git a/Detector/DetCRD/src/Other/CRDBeamPipe_v01_geo.cpp b/Detector/DetCRD/src/Other/CRDBeamPipe_v01_geo.cpp index 8a5c331709de51c0a378ea525b0c2cb7482b9c06..1b1b2a813e716632189dd6b2d54e7df3ade4ad54 100644 --- a/Detector/DetCRD/src/Other/CRDBeamPipe_v01_geo.cpp +++ b/Detector/DetCRD/src/Other/CRDBeamPipe_v01_geo.cpp @@ -2,12 +2,10 @@ // CepC BeamPipe models in DD4hep //-------------------------------------------------------------------- #include "OtherDetectorHelpers.h" -#include "TGeoScaledShape.h" #include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DD4hepUnits.h" #include "DD4hep/DetType.h" -#include "DD4hep/detail/Handle.inl" #include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" @@ -28,8 +26,7 @@ using dd4hep::rec::ConicalSupportData; using dd4hep::rec::SurfaceType; using dd4hep::rec::Vector3D; using dd4hep::rec::VolCylinder; -using dd4hep::rec::VolCylinderImpl; -using dd4hep::rec::VolSurface; +using dd4hep::rec::VolCone; using dd4hep::rec::volSurfaceList; //BeamPipe @@ -62,7 +59,7 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::xml::Component xmlParameter = x_beampipe.child(_Unicode(parameter)); const double crossingAngle = xmlParameter.attr< double >(_Unicode(crossingangle)); std::cout << "Crossing angle = " << crossingAngle << std::endl; - + std::cout << "Section: Zstart Zend RiStart RiEnd size shift type" << std::endl; for(xml_coll_t si( x_beampipe ,Unicode("section")); si; ++si) { xml_comp_t x_section(si); @@ -78,7 +75,7 @@ static Ref_t create_detector(Detector& theDetector, catch(std::runtime_error& e){ rInnerEnd = rInnerStart; } - if(type==CEPC::kWaist || type==CEPC::kFatWaist || type==CEPC::kCrotch || type==CEPC::kCrotchAsymUp || type==CEPC::kCrotchAsymDn){ + if(type==CEPC::kWaist || type==CEPC::kFatWaist || type==CEPC::kCrotch || type==CEPC::kCrotchAsymUp || type==CEPC::kCrotchAsymDn || type==CEPC::kRunway){ try{ size = x_section.attr< double > (_Unicode(size)); } @@ -94,15 +91,15 @@ static Ref_t create_detector(Detector& theDetector, } const std::string volName = "BeamPipe_" + x_section.nameStr(); - - std::cout << "section: " + std::cout << std::setiosflags(std::ios::left) + << std::setw(35) << volName << std::setw(8) << zstart /dd4hep::mm << std::setw(8) << zend /dd4hep::mm << std::setw(8) << rInnerStart /dd4hep::mm << std::setw(8) << rInnerEnd /dd4hep::mm << std::setw(8) << size /dd4hep::mm + << std::setw(8) << shift /dd4hep::mm << std::setw(8) << type - << std::setw(35) << volName << std::endl; const double angle = crossingAngle; @@ -149,7 +146,7 @@ static Ref_t create_detector(Detector& theDetector, catch(std::runtime_error& e){ thicknessEnd = thickness; } - std::cout << " layer: " << std::setw(8) << thickness/dd4hep::mm << std::setw(8) << thicknessEnd/dd4hep::mm << std::setw(15) << material.name() << std::endl; + std::cout << "->layer: " << std::setw(6) << thickness/dd4hep::mm << std::setw(6) << thicknessEnd/dd4hep::mm << std::setw(15) << material.name() << std::endl; char suffix[20]; sprintf(suffix,"_%d",ilayer); @@ -157,13 +154,11 @@ static Ref_t create_detector(Detector& theDetector, if(type==CEPC::kCenter || type==CEPC::kCenterSide){ dd4hep::ConeSegment subLayer(zHalf, radius, radius+thickness, radiusEnd, radiusEnd+thicknessEnd, phi0, dPhi); dd4hep::Volume subLayerLog(volName, subLayer, material); + subLayerLog.setVisAttributes(theDetector, x_layer.visStr()); dd4hep::Transform3D transformer(dd4hep::RotationY(0), dd4hep::Position(0, 0, zCenter)); dd4hep::Transform3D transmirror(dd4hep::RotationY(180*dd4hep::degree), dd4hep::RotateY(dd4hep::Position(0, 0, zCenter), 180*dd4hep::degree)); envelope.placeVolume(subLayerLog, transformer); envelope.placeVolume(subLayerLog, transmirror); - std::cout << "fucd debug: radL = " << material.radLength()/dd4hep::mm << " intL = " << material.intLength()/dd4hep::mm << std::endl; - if(material.radLength()<10000*dd4hep::mm) subLayerLog.setVisAttributes(theDetector, "TubeVis"); - else subLayerLog.setVisAttributes(theDetector, "VacVis"); if(material.radLength()<10000*dd4hep::mm){ double tEff = thickness/material.radLength()*theDetector.material("G4_Be").radLength(); @@ -174,14 +169,30 @@ static Ref_t create_detector(Detector& theDetector, pipeThicknessEnd += tEffEnd; pipeThicknessRel += thickness; pipeThicknessRelEnd += thicknessEnd; + + // only the center pipe and the link are added to surface list for tracking + /* have been added through envelope + if(radius==radiusEnd&&thickness==thicknessEnd){ + Vector3D ocyl(radius+0.5*thickness, 0., 0.); + VolCylinder surf(subLayerLog, SurfaceType(SurfaceType::Helper), 0.5*thickness, 0.5*thickness, ocyl ) ; + volSurfaceList(tube)->push_back(surf); + } + else{ + Vector3D ocon(radius+0.5*thickness, 0., 0.); + Vector3D con_angle(1., 0., atan((radiusEnd+0.5*thicknessEnd-radius-0.5*thickness)/(2.*zHalf)), Vector3D::spherical); + VolCone surf(subLayerLog, SurfaceType(SurfaceType::Helper), 0.5*(thickness+thicknessEnd)/2, 0.5*(thickness+thicknessEnd)/2, con_angle, ocon); + volSurfaceList(tube)->push_back(surf); + } + */ } } - else if(type==CEPC::kLegs){ + else if(type==CEPC::kLegs){ // doubly pipe, same for up and down double clipAngle = 0.5*angle; double lowNorml[3] = { sin(clipAngle), 0, -cos(clipAngle)}; double highNorml[3] = {-sin(clipAngle), 0, cos(clipAngle)}; dd4hep::CutTube subLayer(radius, radius+thickness, zHalf/cos(clipAngle), phi0, dPhi, lowNorml[0], lowNorml[1], lowNorml[2], highNorml[0], highNorml[1], highNorml[2]); dd4hep::Volume subLayerLog(volName, subLayer, material); + subLayerLog.setVisAttributes(theDetector, x_layer.visStr()); dd4hep::Transform3D plusUpTransformer(dd4hep::RotationY(clipAngle), dd4hep::RotateY(dd4hep::Position(0,0,zCenter/cos(clipAngle)), clipAngle)); dd4hep::Transform3D plusDownTransformer(dd4hep::RotationZYX(180*dd4hep::degree, -clipAngle, 0), dd4hep::RotateY(dd4hep::Position(0,0,zCenter/cos(clipAngle)), -clipAngle)); dd4hep::Transform3D minusUpTransformer(dd4hep::RotationY(clipAngle), dd4hep::RotateY(dd4hep::Position(0,0,-zCenter/cos(clipAngle)), clipAngle)); @@ -190,11 +201,8 @@ static Ref_t create_detector(Detector& theDetector, envelope.placeVolume(subLayerLog, plusDownTransformer); envelope.placeVolume(subLayerLog, minusUpTransformer); envelope.placeVolume(subLayerLog, minusDownTransformer); - - if(material.radLength()<10000*dd4hep::mm) subLayerLog.setVisAttributes(theDetector, "TubeVis"); - else subLayerLog.setVisAttributes(theDetector, "VacVis"); } - else if(type==CEPC::kFlareLegUp || type==CEPC::kFlareLegDn){ + else if(type==CEPC::kFlareLegUp || type==CEPC::kFlareLegDn){ // doubly pipe, different for up and down double clipAngle = (type==CEPC::kFlareLegUp)?0.5*angle:-0.5*angle; double rOuter = radius+thickness; double rOuterEnd = radiusEnd+thicknessEnd; @@ -205,14 +213,12 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::ConeSegment wholeSolid(zHalf + clipSize, radius, rOuter, radiusEnd, rOuterEnd, phi0, dPhi); dd4hep::IntersectionSolid layerSolid(wholeSolid, clipSolid, clipTransformer); dd4hep::Volume subLayerLog(volName, layerSolid, material); + subLayerLog.setVisAttributes(theDetector, x_layer.visStr()); envelope.placeVolume(subLayerLog, placementTransformer); envelope.placeVolume(subLayerLog, placementTransmirror); - - if(material.radLength()<10000*dd4hep::mm) subLayerLog.setVisAttributes(theDetector, "TubeVis"); - else subLayerLog.setVisAttributes(theDetector, "VacVis"); } - else if(type==CEPC::kCrotch){ + else if(type==CEPC::kCrotch){ // runway to doubly, center has obstruct between two pipes double beamAngle = 0.5*angle; if(size==0) size = (zstart*tan(beamAngle)+radius)*2; double x1 = 0.5*size - radius; @@ -235,6 +241,7 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::UnionSolid tmp2Solid(tmp1Solid, side1, unionTransformer2); dd4hep::IntersectionSolid shell(tmp2Solid, cut1, sameTransformer); dd4hep::Volume shellLog(volName+"Shell", shell, material); + shellLog.setVisAttributes(theDetector, x_layer.visStr()); envelope.placeVolume(shellLog, dd4hep::Position(0, 0, zCenter)); envelope.placeVolume(shellLog, dd4hep::Transform3D(dd4hep::RotationY(180*dd4hep::degree), dd4hep::Position(0, 0, -zCenter))); @@ -254,12 +261,10 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::UnionSolid tmp7Solid(tmp6Solid, side3, unionTransformer4); dd4hep::IntersectionSolid vacuumPipe(tmp7Solid, cut1, sameTransformer); dd4hep::Volume pipeLog(volName+"Vacuum", vacuumPipe, beamMaterial); + pipeLog.setVisAttributes(theDetector, x_beampipe.visStr()); shellLog.placeVolume(pipeLog, dd4hep::Position(0, 0, 0)); - - shellLog.setVisAttributes(theDetector, "TubeVis"); - pipeLog.setVisAttributes(theDetector, "VacVis"); } - else if(type==CEPC::kCrotchAsymUp || type==CEPC::kCrotchAsymDn){ + else if(type==CEPC::kCrotchAsymUp || type==CEPC::kCrotchAsymDn){ // runway to doubly, center has obstruct between two pipes, and different between up and down double beamAngle = 0.5*angle; double xC2 = (shift==0)?zend*tan(beamAngle):shift; if(radiusEnd==0) radiusEnd = radius; @@ -301,11 +306,8 @@ static Ref_t create_detector(Detector& theDetector, double thetaCut1 = atan((0.5*(xC2+rOuterEnd)-0.5*rOuter)/(2*zHalf)); double xcenterCut1 = 0.5*(0.5*(xC2+rOuterEnd)+0.5*rOuter); dd4hep::Trap cut1(zHalf, thetaCut1, 0, yMax, 0.5*rOuter, 0.5*rOuter, 0, rOuterEnd, 0.5*(xC2+rOuterEnd), 0.5*(xC2+rOuterEnd), 0); - TGeoCone* pCone1 = new TGeoCone(pzTopCut, 0, a1, 0, a2); - //double factor = - TGeoScale* pScale1 = new TGeoScale(1, b1/a1, 1); - TGeoScaledShape* pScaledShape1 = new TGeoScaledShape(pCone1, pScale1); - dd4hep::Solid_type<TGeoScaledShape> side1(pScaledShape1); + dd4hep::Cone cone1(pzTopCut, 0, a1, 0, a2); + dd4hep::Scale side1(cone1, 1, b1/a1, 1); double xshift = 0.5*(xMaxEnd-a2*cos(rotateAngle)-rOuter+a1*cos(bottomAngle-edge2ToXAngle)); double zshift = 0.5*(a2-a1)*sin(rotateAngle); @@ -314,6 +316,7 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::UnionSolid tmp1Solid(body, side1, unionTransformer1); dd4hep::IntersectionSolid shell(tmp1Solid, cut1, cutTransformer1); dd4hep::Volume shellLog(volName, shell, material); + shellLog.setVisAttributes(theDetector, x_layer.visStr()); if(type==CEPC::kCrotchAsymUp){ envelope.placeVolume(shellLog, dd4hep::Position(0, 0, zCenter)); envelope.placeVolume(shellLog, dd4hep::Transform3D(dd4hep::RotationX(180*dd4hep::degree), dd4hep::Position(0, 0, -zCenter))); @@ -353,21 +356,17 @@ static Ref_t create_detector(Detector& theDetector, double thetaCut2 = atan((xC2-0.5*radius)/(2*zHalf)); double xcenterCut2 = 0.5*radius+0.5*(xC2-0.5*radius); dd4hep::Trap cut2(zHalf, thetaCut2, 0, yMax-thickness, 0.5*radius, 0.5*radius, 0, radiusEnd, radiusEnd, radiusEnd, 0); - TGeoCone* pCone2 = new TGeoCone(pzTopCutHole, 0, a1Hole, 0, a2Hole); - TGeoScale* pScale2 = new TGeoScale(1, b1Hole/a1Hole, 1); - TGeoScaledShape* pScaledShape2 = new TGeoScaledShape(pCone2, pScale2); - dd4hep::Solid_type<TGeoScaledShape> side2(pScaledShape2); + dd4hep::Cone cone2(pzTopCutHole, 0, a1Hole, 0, a2Hole); + dd4hep::Scale side2(cone2, 1, b1Hole/a1Hole, 1); double xshiftHole = 0.5*(xMaxEnd-thicknessEnd-a2Hole*cos(rotate)-radius+a1Hole*cos(bottom-edge2ToX)); double zshiftHole = 0.5*(a2Hole-a1Hole)*sin(rotate); dd4hep::Transform3D cutTransformer2(dd4hep::RotationY(rotate), dd4hep::Position(xshiftHole-xcenterCut2, 0, zshiftHole)); dd4hep::IntersectionSolid vacuumPipe(cut2, side2, cutTransformer2); dd4hep::Volume pipeLog(volName, vacuumPipe, beamMaterial); + pipeLog.setVisAttributes(theDetector, x_beampipe.visStr()); shellLog.placeVolume(pipeLog, dd4hep::Position(xcenterCut2, 0, 0)); - - shellLog.setVisAttributes(theDetector, "TubeVis"); - pipeLog.setVisAttributes(theDetector, "VacVis"); } - else if(type==CEPC::kWaist){ + else if(type==CEPC::kWaist){ // expanded single pipe from circle to runway (each turn: 180 degree) double beamAngle = 0.5*angle; if(radiusEnd==0) radiusEnd = radius; if(size==0) size = (zend*tan(beamAngle)+radiusEnd)*2; @@ -379,37 +378,42 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::Trd2 body1(0, xC2, rOuter, rOuterEnd, zHalf); dd4hep::Trd2 cut(rOuter, rMaxEnd, rOuter, rOuterEnd, zHalf); - double expandAngle = atan(xC2/(2*zHalf)); - double edge1ToZAngle = atan((rMaxEnd-rOuter)/(2*zHalf)); - double edge2ToZAngle = atan((xC2-rOuterEnd+rOuter)/(2*zHalf)); - double edge2ToXAngle = 90*dd4hep::degree - edge2ToZAngle; - double bottomAngle = 0.5*(180*dd4hep::degree-(edge2ToZAngle-edge1ToZAngle)); - double rotateAngle = 0.5*(edge1ToZAngle+edge2ToZAngle); - double edge1ToCAngle = asin(sin(90*dd4hep::degree+edge1ToZAngle)/(xC2/sin(expandAngle))*(rOuter-rOuterEnd)); - double CToEConeAxisAngle = edge1ToCAngle-0.5*(edge2ToZAngle-edge1ToZAngle); - if(fabs(rotateAngle-(expandAngle-CToEConeAxisAngle))>1e-12){ + double pzTopCut = zHalf; + double expandAngle = atan(xC2/(2*zHalf)); + double edge1ToZAngle = atan((rMaxEnd-rOuter)/(2*zHalf)); + double edge2ToZAngle = atan((xC2-rOuterEnd+rOuter)/(2*zHalf)); + double edge2ToXAngle = 90*dd4hep::degree - edge2ToZAngle; + double bottomAngle = 0.5*(180*dd4hep::degree-(edge2ToZAngle-edge1ToZAngle)); + double rotateAngle = 0.5*(edge1ToZAngle+edge2ToZAngle); + double edge1ToCAngle = asin(sin(90*dd4hep::degree+edge1ToZAngle)/(xC2/sin(expandAngle))*(rOuter-rOuterEnd)); + double CToEConeAxisAngle = edge1ToCAngle-0.5*(edge2ToZAngle-edge1ToZAngle); + std::cout << expandAngle/dd4hep::degree << " " << edge1ToZAngle/dd4hep::degree << " " << edge2ToZAngle/dd4hep::degree << " " << bottomAngle/dd4hep::degree << " " << rotateAngle/dd4hep::degree << " " << edge1ToCAngle/dd4hep::degree << " " << CToEConeAxisAngle/dd4hep::degree << std::endl; + if(fabs(rotateAngle-(expandAngle-CToEConeAxisAngle))>1e-12){ std::cout << "Warning! rotate angle was not calculated rightly. Please check input parameters whether satisfy the Waist case." << std::endl; - } + } double a1 = rOuter/sin(bottomAngle)*sin(90*dd4hep::degree-edge1ToZAngle); - double a2 = rOuterEnd/sin(180*dd4hep::degree-bottomAngle)*sin(90*dd4hep::degree-edge2ToZAngle); - double zC1 = rOuter/sin(edge1ToCAngle)*sin(90*dd4hep::degree+edge1ToZAngle)*cos(CToEConeAxisAngle); - double zC2 = rOuterEnd/rOuter*zC1; - double zBottom = a1*tan(bottomAngle); - double aC1 = a1/zBottom*zC1; - double aC2 = a1/zBottom*zC2; - double xC1InECone = zC1*tan(CToEConeAxisAngle); - double xC2InECone = zC2*tan(CToEConeAxisAngle); - double bC1 = sqrt(rOuter*rOuter/(1-xC1InECone*xC1InECone/aC1/aC1)); - double bC2 = sqrt(rOuterEnd*rOuterEnd/(1-xC2InECone*xC2InECone/aC2/aC2)); - double b1 = bC1/zC1*zBottom; - if(fabs(bC1/zC1-bC2/zC2)>1e-12){ - std::cout << "Warning! bC1/zC1 not equal to bC2/zC2." << std::endl; - } - double pzTopCut = 0.5*(a1-a2)*tan(bottomAngle); - TGeoCone* pCone1 = new TGeoCone(pzTopCut, 0, a1, 0, a2); - TGeoScale* pScale1 = new TGeoScale(1, b1/a1, 1); - TGeoScaledShape* pScaledShape1 = new TGeoScaledShape(pCone1,pScale1); - dd4hep::Solid_type<TGeoScaledShape> side1(pScaledShape1); + double a2 = rOuterEnd/sin(180*dd4hep::degree-bottomAngle)*sin(90*dd4hep::degree-edge2ToZAngle); + double b1 = rOuter; + if(rOuter!=rOuterEnd){ + double zC1 = rOuter/sin(edge1ToCAngle)*sin(90*dd4hep::degree+edge1ToZAngle)*cos(CToEConeAxisAngle); + double zC2 = rOuterEnd/rOuter*zC1; + double zBottom = a1*tan(bottomAngle); + double aC1 = a1/zBottom*zC1; + double aC2 = a1/zBottom*zC2; + double xC1InECone = zC1*tan(CToEConeAxisAngle); + double xC2InECone = zC2*tan(CToEConeAxisAngle); + double bC1 = sqrt(rOuter*rOuter/(1-xC1InECone*xC1InECone/aC1/aC1)); + double bC2 = sqrt(rOuterEnd*rOuterEnd/(1-xC2InECone*xC2InECone/aC2/aC2)); + b1 = bC1/zC1*zBottom; + std::cout << a1 << " " << a2 << " " << zC1 << " " << zC2 << std::endl; + if(fabs(bC1/zC1-bC2/zC2)>1e-12){ + std::cout << "Warning! bC1/zC1 not equal to bC2/zC2." << std::endl; + } + std::cout << "b1/a1=" << b1/a1 << std::endl; + pzTopCut = 0.5*(a1-a2)*tan(bottomAngle); + } + dd4hep::Cone cone1(pzTopCut, 0, a1, 0, a2); + dd4hep::Scale side1(cone1, 1, b1/a1, 1); double xshift = 0.5*(rMaxEnd-a2*cos(rotateAngle)-rOuter+a1*cos(bottomAngle-edge2ToXAngle)); double zshift = 0.5*(a2-a1)*sin(rotateAngle); @@ -420,9 +424,11 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::UnionSolid tmp2Solid(tmp1Solid, side1, unionTransformer2); dd4hep::IntersectionSolid shell(tmp2Solid, cut, sameTransformer); dd4hep::Volume shellLog(volName+"Shell", shell, material); + shellLog.setVisAttributes(theDetector, x_layer.visStr()); envelope.placeVolume(shellLog, dd4hep::Position(0, 0, zCenter)); envelope.placeVolume(shellLog, dd4hep::Transform3D(dd4hep::RotationY(180*dd4hep::degree), dd4hep::Position(0, 0, -zCenter))); - + + double pzTopCutHole = zHalf; double edge1ToZ = atan((0.5*size-radius)/(2*zHalf)); double edge2ToZ = atan((xC2-radiusEnd+radius)/(2*zHalf)); double edge2ToX = 90*dd4hep::degree - edge2ToZ; @@ -435,26 +441,27 @@ static Ref_t create_detector(Detector& theDetector, } double a1Hole = radius/sin(bottom)*sin(90*dd4hep::degree-edge1ToZ); double a2Hole = radiusEnd/sin(180*dd4hep::degree-bottom)*sin(90*dd4hep::degree-edge2ToZ); - double zC1Hole = radius/sin(edge1ToC)*sin(90*dd4hep::degree+edge1ToZ)*cos(CToEConeAxis); - double zC2Hole = radiusEnd/radius*zC1Hole; - double zBottomHole = a1Hole*tan(bottom); - double aC1Hole = a1Hole/zBottomHole*zC1Hole; - double aC2Hole = a1Hole/zBottomHole*zC2Hole; - double xC1InEConeHole = zC1Hole*tan(CToEConeAxis); - double xC2InEConeHole = zC2Hole*tan(CToEConeAxis); - double bC1Hole = sqrt(radius*radius/(1-xC1InEConeHole*xC1InEConeHole/aC1Hole/aC1Hole)); - double bC2Hole = sqrt(radiusEnd*radiusEnd/(1-xC2InEConeHole*xC2InEConeHole/aC2Hole/aC2Hole)); - double b1Hole = bC1Hole/zC1Hole*zBottomHole; - if(fabs(bC1Hole/zC1Hole-bC2Hole/zC2Hole)>1e-12){ - std::cout << "Warning! bC1/zC1 not equal to bC2/zC2 for Hole." << std::endl; - } - double pzTopCutHole = 0.5*(a1Hole-a2Hole)*tan(bottom); + double b1Hole = radius; + if(radius!=radiusEnd){ + double zC1Hole = radius/sin(edge1ToC)*sin(90*dd4hep::degree+edge1ToZ)*cos(CToEConeAxis); + double zC2Hole = radiusEnd/radius*zC1Hole; + double zBottomHole = a1Hole*tan(bottom); + double aC1Hole = a1Hole/zBottomHole*zC1Hole; + double aC2Hole = a1Hole/zBottomHole*zC2Hole; + double xC1InEConeHole = zC1Hole*tan(CToEConeAxis); + double xC2InEConeHole = zC2Hole*tan(CToEConeAxis); + double bC1Hole = sqrt(radius*radius/(1-xC1InEConeHole*xC1InEConeHole/aC1Hole/aC1Hole)); + double bC2Hole = sqrt(radiusEnd*radiusEnd/(1-xC2InEConeHole*xC2InEConeHole/aC2Hole/aC2Hole)); + b1Hole = bC1Hole/zC1Hole*zBottomHole; + if(fabs(bC1Hole/zC1Hole-bC2Hole/zC2Hole)>1e-12){ + std::cout << "Warning! bC1/zC1 not equal to bC2/zC2 for Hole." << std::endl; + } + pzTopCutHole = 0.5*(a1Hole-a2Hole)*tan(bottom); + } dd4hep::Trd2 body2(0, xC2, radius, radiusEnd, zHalf); dd4hep::Trd2 cut2(radius, 0.5*size, radius, radiusEnd, zHalf); - TGeoCone* pCone2 = new TGeoCone(pzTopCutHole, 0, a1Hole, 0, a2Hole); - TGeoScale* pScale2 = new TGeoScale(1, b1Hole/a1Hole, 1); - TGeoScaledShape* pScaledShape2 = new TGeoScaledShape(pCone2,pScale2); - dd4hep::Solid_type<TGeoScaledShape> side2(pScaledShape2); + dd4hep::Cone cone2(pzTopCutHole, 0, a1Hole, 0, a2Hole); + dd4hep::Scale side2(cone2, 1, b1Hole/a1Hole, 1); double xshiftHole = 0.5*(0.5*size-a2Hole*cos(rotate)-radius+a1Hole*cos(bottom-edge2ToX)); double zshiftHole = 0.5*(a2Hole-a1Hole)*sin(rotate); @@ -464,12 +471,10 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::UnionSolid tmp4Solid(tmp3Solid, side2, unionTransformer4); dd4hep::IntersectionSolid vacuumPipe(tmp4Solid, cut, sameTransformer); dd4hep::Volume pipeLog(volName+"Vacuum", vacuumPipe, beamMaterial); + pipeLog.setVisAttributes(theDetector, x_beampipe.visStr()); shellLog.placeVolume(pipeLog, dd4hep::Position(0, 0, 0)); - - shellLog.setVisAttributes(theDetector, "TubeVis"); - pipeLog.setVisAttributes(theDetector, "VacVis"); } - else if(type == CEPC::kFatWaist){ + else if(type == CEPC::kFatWaist){ // expanded single pipe from circle to cuted circle, runway but not 180 degree double beamAngle = 0.5*angle; if(radiusEnd==0) radiusEnd = radius; if(size==0) size = (zend*tan(beamAngle)+radiusEnd)*2; @@ -483,6 +488,7 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::ConeSegment cone1(zHalf, 0, rOuter, 0, rOuterEnd, phi0, dPhi); dd4hep::IntersectionSolid shell(cone1, body1, sameTransformer); dd4hep::Volume shellLog(volName, shell, material); + shellLog.setVisAttributes(theDetector, x_layer.visStr()); envelope.placeVolume(shellLog, dd4hep::Position(0, 0, zCenter)); envelope.placeVolume(shellLog, dd4hep::Transform3D(dd4hep::RotationY(180*dd4hep::degree), dd4hep::Position(0, 0, -zCenter))); @@ -490,10 +496,47 @@ static Ref_t create_detector(Detector& theDetector, dd4hep::ConeSegment cone2(zHalf, 0, radius, 0, radius, phi0, dPhi); dd4hep::IntersectionSolid vacuumPipe(cone2, body2, sameTransformer); dd4hep::Volume pipeLog(volName, vacuumPipe, beamMaterial); + pipeLog.setVisAttributes(theDetector, x_beampipe.visStr()); shellLog.placeVolume(pipeLog, dd4hep::Position(0, 0, 0)); + } + else if(type == CEPC::kRunway){ // runway to runway (180 degree), same radius + double sizeEnd = size + shift; + if(size==0){ + size = (zstart*tan(0.5*angle)+radius)*2; + sizeEnd = (zend*tan(0.5*angle)+radius)*2; + } + double x1 = 0.5*size - radius; + double x2 = 0.5*sizeEnd - radius; + double y1 = radius + thickness; + double y2 = y1; + double expandAngle = atan((x2-x1)/zHalf/2); + double zSide = 2*zHalf/cos(expandAngle)+y1*tan(expandAngle)+y2*tan(expandAngle); + double xshift = 0.5*(x1+x2); - shellLog.setVisAttributes(theDetector, "TubeVis"); - pipeLog.setVisAttributes(theDetector, "VacVis"); + dd4hep::Trd2 body1(x1, x2, y1, y2, zHalf); + dd4hep::Trd2 cut1(x1+y1, x2+y2, y1, y2, zHalf); + dd4hep::EllipticalTube side1(y1*cos(expandAngle), y1, 0.5*zSide); + dd4hep::Transform3D unionTransformer1(dd4hep::RotationY(expandAngle), dd4hep::Position(xshift, 0, 0)); + dd4hep::Transform3D unionTransformer2(dd4hep::RotationY(-expandAngle), dd4hep::Position(-xshift, 0, 0)); + dd4hep::Transform3D sameTransformer(dd4hep::RotationY(0), dd4hep::Position(0, 0, 0)); + dd4hep::UnionSolid tmp1Solid(body1, side1, unionTransformer1); + dd4hep::UnionSolid tmp2Solid(tmp1Solid, side1, unionTransformer2); + dd4hep::IntersectionSolid shell(tmp2Solid, cut1, sameTransformer); + dd4hep::Volume shellLog(volName+"Shell", shell, material); + shellLog.setVisAttributes(theDetector, x_layer.visStr()); + envelope.placeVolume(shellLog, dd4hep::Position(0, 0, zCenter)); + envelope.placeVolume(shellLog, dd4hep::Transform3D(dd4hep::RotationY(180*dd4hep::degree), dd4hep::Position(0, 0, -zCenter))); + + double yHole = y1-thickness; + dd4hep::Trd2 body2(x1, x2, yHole, yHole, zHalf); + dd4hep::Trd2 cut2(x1+y1, x2+y2, yHole, yHole, zHalf); + dd4hep::EllipticalTube side2(yHole*cos(expandAngle), yHole, zSide); + dd4hep::UnionSolid tmp3Solid(body2, side2, unionTransformer1); + dd4hep::UnionSolid tmp4Solid(tmp3Solid, side2, unionTransformer2); + dd4hep::IntersectionSolid vacuumPipe(tmp4Solid, cut2, sameTransformer); + dd4hep::Volume pipeLog(volName+"Vacuum", vacuumPipe, beamMaterial); + pipeLog.setVisAttributes(theDetector, x_beampipe.visStr()); + shellLog.placeVolume(pipeLog, dd4hep::Position(0, 0, 0)); } radius += thickness; radiusEnd += thicknessEnd; @@ -532,7 +575,7 @@ static Ref_t create_detector(Detector& theDetector, tube.addExtension< ConicalSupportData >( beampipeData ) ; //-------------------------------------- - tube.setVisAttributes( theDetector, x_beampipe.visStr(), envelope ); + tube.setVisAttributes( theDetector, "SeeThrough", envelope ); //debug std::cout << "============ConicalSupportData============" << std::endl; @@ -545,6 +588,4 @@ static Ref_t create_detector(Detector& theDetector, return tube; } -DECLARE_DETELEMENT(DD4hep_CRDBeamPipe_v01, create_detector) - -DD4HEP_INSTANTIATE_SHAPE_HANDLE(TGeoScaledShape); +DECLARE_DETELEMENT(CRDBeamPipe_v01, create_detector) diff --git a/Detector/DetCRD/src/Other/OtherDetectorHelpers.h b/Detector/DetCRD/src/Other/OtherDetectorHelpers.h index 7e7569af39dce2667afad5ff75dcf58813a7f5e9..86efbfccf32d5c1dc41fc4618821a7dbdfa12772 100644 --- a/Detector/DetCRD/src/Other/OtherDetectorHelpers.h +++ b/Detector/DetCRD/src/Other/OtherDetectorHelpers.h @@ -16,7 +16,8 @@ namespace CEPC { kCrotchAsymDn = 6, kLegs = 7, kFlareLegUp = 8, - kFlareLegDn = 9 + kFlareLegDn = 9, + kRunway = 10 } ECrossType; inline ECrossType getCrossType( std::string const & type) { @@ -30,7 +31,8 @@ namespace CEPC { CrossTypes["CrotchAsymDn"] = CEPC::kCrotchAsymDn; CrossTypes["Legs"] = CEPC::kLegs; CrossTypes["FlareLegUp"] = CEPC::kFlareLegUp; - CrossTypes["FlareLegDn"] = CEPC::kFlareLegDn; + CrossTypes["FlareLegDn"] = CEPC::kFlareLegDn; + CrossTypes["Runway"] = CEPC::kRunway; std::map < std::string, CEPC::ECrossType>::const_iterator it = CrossTypes.find(type); if ( it == CrossTypes.end() ) { diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1496aff196289e0d5700f48255054a38bcc5f3f9 --- /dev/null +++ b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp @@ -0,0 +1,402 @@ +//==================================================================== +// cepcvxdgeo - CEPC vertex detector models in DD4hep +//-------------------------------------------------------------------- +// Hao Zeng, IHEP +// email: zenghao@ihep.ac.cn +// $Id$ +//==================================================================== +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/DetType.h" +#include "DDRec/Surface.h" +#include "DDRec/DetectorData.h" +#include "XML/Utilities.h" +#include <cmath> + +using namespace std; + +using dd4hep::Box; +using dd4hep::DetElement; +using dd4hep::Material; +using dd4hep::Position; +using dd4hep::RotationY; +using dd4hep::RotationZYX; +using dd4hep::Transform3D; +using dd4hep::Rotation3D; +using dd4hep::Volume; +using dd4hep::_toString; +using dd4hep::rec::volSurfaceList; +using dd4hep::rec::ZPlanarData; +using dd4hep::mm; + +/** helper struct */ +struct VXD_Layer { + int n_ladders; + int n_sensors_per_ladder; + double sensor_length; + double half_z; + double sensitive_inner_radius ; + double support_inner_radius ; + double ladder_width ; + double ladder_dphi ; +}; + +//std::vector<VXD_Layer> _VXD_Layers; + +// /** helper struct */ +// struct extended_reconstruction_parameters { +// double sensor_length_mm; +// double strip_width_mm; +// double strip_length_mm; +// double strip_pitch_mm; +// double strip_angle_deg; +// }; + +//extended_reconstruction_parameters _e_r_p; + + +/** Construction of the VXD detector, ported from Mokka driver SIT_Simple_Pixel.cc + * + * Mokka History: + * Feb 7th 2011, Steve Aplin - original version + * F.Gaede, DESY, Jan 2014 - dd4hep SIT pixel + + * @author Hao Zeng, IHEP, July 2021 + */ +static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4hep::SensitiveDetector sens) { + + xml_det_t x_det = e; + Material air = theDetector.air(); + int det_id = x_det.id(); + string name = x_det.nameStr(); + DetElement vxd(name, det_id); + + Volume envelope = dd4hep::xml::createPlacedEnvelope(theDetector, e, vxd); + dd4hep::xml::setDetectorTypeFlag(e, vxd) ; + if(theDetector.buildType()==dd4hep::BUILD_ENVELOPE) return vxd; + envelope.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + sens.setType("tracker"); + std::cout << " ** building SiTrackerSkewBarrel_v01 ..." << std::endl ; + + dd4hep::rec::ZPlanarData* zPlanarData = new dd4hep::rec::ZPlanarData; + + // fetch the global parameters + + //fetch the display parameters + xml_comp_t x_display(x_det.child(_Unicode(display))); + std::string ladderVis = x_display.attr<string>(_Unicode(ladder)); + std::string supportVis = x_display.attr<string>(_Unicode(support)); + std::string flexVis = x_display.attr<string>(_Unicode(flex)); + std::string sensEnvVis = x_display.attr<string>(_Unicode(sens_env)); + std::string sensVis = x_display.attr<string>(_Unicode(sens)); + std::string deadsensVis = x_display.attr<string>(_Unicode(deadsensor)); + std::string deadwireVis = x_display.attr<string>(_Unicode(deadwire)); + + + for(xml_coll_t layer_i(x_det,_U(layer)); layer_i; ++layer_i){ + xml_comp_t x_layer(layer_i); + + dd4hep::PlacedVolume pv; + int layer_id = x_layer.attr<int>(_Unicode(layer_id)); + + std::cout << "layer_id: " << layer_id << endl; + + double sensitive_radius = x_layer.attr<double>(_Unicode(ladder_radius)); + int n_sensors_per_ladder = x_layer.attr<int>(_Unicode(n_sensors_per_side)); + int n_ladders = x_layer.attr<int>(_Unicode(n_ladders)) ; + double ladder_offset = x_layer.attr<double>(_Unicode(ladder_offset)); + double ladder_radius = sqrt(ladder_offset*ladder_offset + sensitive_radius*sensitive_radius); + double ladder_phi0 = -atan(ladder_offset/sensitive_radius); + std::cout << "ladder_radius: " << ladder_radius/mm <<" mm" << endl; + + std::cout << "sensitive_radius: " << sensitive_radius/mm << " mm" << endl; + std::cout << "n_sensors_per_ladder: " << n_sensors_per_ladder << endl; + + std::string layerName = dd4hep::_toString( layer_id , "layer_%d" ); + dd4hep::Assembly layer_assembly( layerName ) ; + pv = envelope.placeVolume( layer_assembly ) ; + dd4hep::DetElement layerDE( vxd , layerName , x_det.id() ); + layerDE.setPlacement( pv ) ; + + const double ladder_dphi = ( dd4hep::twopi / n_ladders ) ; + std::cout << "ladder_dphi: " << ladder_dphi << endl; + + //fetch the ladder parameters + xml_comp_t x_ladder(x_layer.child(_Unicode(ladder))); + // db = XMLHandlerDB(x_ladder); + + //fetch the ladder support parameters + xml_comp_t x_ladder_support(x_ladder.child(_Unicode(ladderSupport))); + double support_length = x_ladder_support.attr<double>(_Unicode(length)); + double support_thickness = x_ladder_support.attr<double>(_Unicode(thickness)); + double support_height = x_ladder_support.attr<double>(_Unicode(height)); + double support_width = x_ladder_support.attr<double>(_Unicode(width)); + Material support_mat; + if(x_ladder_support.hasAttr(_Unicode(mat))) + { + support_mat = theDetector.material(x_ladder_support.attr<string>(_Unicode(mat))); + } + else + { + support_mat = theDetector.material(x_ladder_support.materialStr()); + } + std::cout << "support_length: " << support_length/mm << " mm" << endl; + std::cout << "support_thickness: " << support_thickness/mm << " mm" << endl; + std::cout << "support_width: " << support_width/mm << " mm" << endl; + + //fetch the flex parameters + double flex_thickness(0); + double flex_width(0); + double flex_length(0); + xml_comp_t x_flex(x_ladder.child(_Unicode(flex))); + for(xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i){ + xml_comp_t x_flex_slice(flex_i); + double x_flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness)); + double x_flex_slice_width = x_flex_slice.attr<double>(_Unicode(width)); + double x_flex_slice_length = x_flex_slice.attr<double>(_Unicode(length)); + flex_thickness += x_flex_slice_thickness; + if (x_flex_slice_width > flex_width) flex_width = x_flex_slice_width; + if (x_flex_slice_length > flex_length) flex_length = x_flex_slice_length; + std::cout << "x_flex_slice_thickness: " << x_flex_slice_thickness/mm << " mm" << endl; + } + std::cout << "flex_thickness: " << flex_thickness/mm << " mm" << endl; + std::cout << "flex_width: " << flex_width/mm << " mm" << endl; + std::cout << "flex_length: " << flex_length/mm << " mm" << endl; + + //fetch the sensor parameters + xml_comp_t x_sensor(x_ladder.child(_Unicode(sensor))); + int n_sensors_per_side = x_sensor.attr<int>(_Unicode(n_sensors)); + double dead_gap = x_sensor.attr<double>(_Unicode(gap)); + double sensor_thickness = x_sensor.attr<double>(_Unicode(thickness)); + double sensor_active_len = x_sensor.attr<double>(_Unicode(active_length)); + double sensor_active_width = x_sensor.attr<double>(_Unicode(active_width)); + double sensor_dead_width = x_sensor.attr<double>(_Unicode(dead_width)); + double sensor_deadwire_length = x_sensor.attr<double>(_Unicode(deadwire_length)); + double sensor_deadwire_width = x_sensor.attr<double>(_Unicode(deadwire_width)); + double sensor_deadwire_thickness = x_sensor.attr<double>(_Unicode(deadwire_thickness)); + Material sensor_mat = theDetector.material(x_sensor.attr<string>(_Unicode(sensor_mat))); + Material sensor_deadwire_mat = theDetector.material(x_sensor.attr<string>(_Unicode(deadwire_mat))); + + std::cout << "n_sensors_per_side: " << n_sensors_per_side << endl; + std::cout << "dead_gap: " << dead_gap/mm << " mm" << endl; + std::cout << "sensor_thickness: " << sensor_thickness/mm << " mm" << endl; + std::cout << "sensor_active_len: " << sensor_active_len/mm << " mm" << endl; + std::cout << "sensor_active_width: " << sensor_active_width/mm << " mm" << endl; + std::cout << "sensor_dead_width: " << sensor_dead_width/mm << " mm" << endl; + + //create ladder logical volume + Box LadderSolid((support_height+2*sensor_thickness+2*flex_thickness)/2.0, + support_width / 2.0, support_length / 2.0); + Volume LadderLogical(name + dd4hep::_toString( layer_id, "_LadderLogical_%02d"), + LadderSolid, air); + // create flex envelope logical volume + Box FlexEnvelopeSolid(flex_thickness / 2.0, flex_width / 2.0, flex_length / 2.0); + Volume FlexEnvelopeLogical(name + dd4hep::_toString( layer_id, "_FlexEnvelopeLogical_%02d"), FlexEnvelopeSolid, air); + FlexEnvelopeLogical.setVisAttributes(theDetector.visAttributes("SeeThrough")); + //vxd.setVisAttributes(theDetector, flexVis, FlexEnvelopeLogical); + + //create the flex layers inside the flex envelope + double flex_start_height(-flex_thickness/2.); + int index = 0; + for(xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i){ + xml_comp_t x_flex_slice(flex_i); + double x_flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness)); + double x_flex_slice_width = x_flex_slice.attr<double>(_Unicode(width)); + double x_flex_slice_length = x_flex_slice.attr<double>(_Unicode(length)); + Material x_flex_slice_mat; + if(x_flex_slice.hasAttr(_Unicode(mat))) + { + x_flex_slice_mat = theDetector.material(x_flex_slice.attr<string>(_Unicode(mat))); + } + else + { + x_flex_slice_mat = theDetector.material(x_flex_slice.materialStr()); + } + // Material x_flex_slice_mat = theDetector.material(x_flex_slice.attr<string>(_Unicode(mat))); + Box FlexLayerSolid(x_flex_slice_thickness/2.0, x_flex_slice_width/2.0, x_flex_slice_length/2.0); + Volume FlexLayerLogical(name + dd4hep::_toString( layer_id, "_FlexLayerLogical_%02d") + dd4hep::_toString( index, "index_%02d"), FlexLayerSolid, x_flex_slice_mat); + FlexLayerLogical.setVisAttributes(theDetector.visAttributes(flexVis)); + double flex_slice_height = flex_start_height + x_flex_slice_thickness/2.; + pv = FlexEnvelopeLogical.placeVolume(FlexLayerLogical, Position(flex_slice_height, 0., 0.)); + std::cout << "flex thickness = " << x_flex_slice_thickness << std::endl; + std::cout << "flex width = " << x_flex_slice_width << std::endl; + std::cout << "flex length = " << x_flex_slice_length << std::endl; + // std::cout << "flex material: " << x_flex_slice_mat << std::endl; + flex_start_height += x_flex_slice_thickness; + index++; + } + + //place the flex envelope inside the ladder envelope + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, Position((support_height + flex_thickness) / 2.0, 0., 0.)); //top side + //define the transformation3D(only need a combination of translation and rotation) + Transform3D tran_mirro(RotationZYX(0., dd4hep::twopi/2.0, 0.), Position(-(support_height + flex_thickness) / 2.0, 0., 0.)); + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, tran_mirro); //bottom side + + //create sensor envelope logical volume + Box SensorTopEnvelopeSolid(sensor_thickness / 2.0, support_width / 2.0, support_length / 2.0); + Volume SensorTopEnvelopeLogical(name + dd4hep::_toString( layer_id, "_SensorEnvelopeLogical_%02d"), SensorTopEnvelopeSolid, air); + Box SensorBottomEnvelopeSolid(sensor_thickness / 2.0, support_width / 2.0, support_length / 2.0); + Volume SensorBottomEnvelopeLogical(name + dd4hep::_toString( layer_id, "_SensorEnvelopeLogical_%02d"), SensorBottomEnvelopeSolid, air); + SensorTopEnvelopeLogical.setVisAttributes(theDetector.visAttributes(sensEnvVis)); + + //create sensor logical volume + Box SensorSolid(sensor_thickness / 2.0, sensor_active_width / 2.0, sensor_active_len / 2.0); + Volume SensorLogical(name + dd4hep::_toString( layer_id, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + SensorLogical.setSensitiveDetector(sens); + //vxd.setVisAttributes(theDetector, deadsensVis, SensorDeadLogical); + SensorLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + + //create dead sensor logical volume + Box SensorDeadSolid(sensor_thickness / 2.0, sensor_dead_width / 2.0, sensor_active_len / 2.0); + Volume SensorDeadLogical(name + dd4hep::_toString( layer_id, "_SensorDeadLogical_%02d"), SensorDeadSolid, sensor_mat); + SensorDeadLogical.setVisAttributes(theDetector.visAttributes(deadsensVis)); + + //create dead wire logical volume + Box SensorDeadWireSolid(sensor_deadwire_thickness / 2.0, sensor_deadwire_width / 2.0, sensor_deadwire_length / 2.0); + Volume SensorDeadWireLogical(name + dd4hep::_toString( layer_id, "_SensorDeadWireLogical_%02d"), SensorDeadWireSolid, sensor_deadwire_mat); + SensorDeadWireLogical.setVisAttributes(theDetector.visAttributes(deadwireVis)); + + //place the dead wire in the sensor envelope + // pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (sensor_active_width-support_width/2.0) + sensor_dead_width/2.0 + sensor_deadwire_width/2.0, 0.0)); + // pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (sensor_active_width-support_width/2.0) + sensor_dead_width/2.0 + sensor_deadwire_width/2.0, 0.0)); + pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (-support_width/2.0) + (sensor_deadwire_width/2.0), 0.0)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (-support_width/2.0) + (sensor_deadwire_width/2.0), 0.0)); + + // place the active sensor and dead sensor inside the sensor envelope + std::vector<dd4hep::PlacedVolume> TopSensor_pv; + std::vector<dd4hep::PlacedVolume> BottomSensor_pv; + for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ + double sensor_total_z = n_sensors_per_side*sensor_active_len + dead_gap*(n_sensors_per_side-1); + double xpos = 0.0; + double ypos_active = (support_width/2.0) - (sensor_active_width/2.0); + double ypos_dead = (-support_width/2.0) + sensor_deadwire_width + (sensor_dead_width/2.0); + double zpos = -sensor_total_z/2.0 + sensor_active_len/2.0 + isensor*(sensor_active_len + dead_gap); + pv = SensorTopEnvelopeLogical.placeVolume(SensorLogical, Position(xpos,ypos_active,zpos)); + //pv.addPhysVolID("topsensor", isensor ) ; + pv.addPhysVolID("layer", layer_id*2+1).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ; + TopSensor_pv.push_back(pv); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorLogical, Position(xpos,ypos_active,zpos)); + //pv.addPhysVolID("bottomsensor", isensor ) ; + pv.addPhysVolID("layer", layer_id*2 ).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ; + BottomSensor_pv.push_back(pv); + pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos)); + + } + //place the sensor envelope inside the ladder envelope + pv = LadderLogical.placeVolume(SensorTopEnvelopeLogical, + Position(support_height/2.0 + flex_thickness + sensor_thickness/2.0, 0., 0.));//top-side sensors + Transform3D tran_sen(RotationZYX(0., dd4hep::twopi/2.0, 0.), Position(-(support_height/2.0 + flex_thickness + sensor_thickness/2.0), 0., 0.)); + pv = LadderLogical.placeVolume(SensorBottomEnvelopeLogical,tran_sen);//bottom-side sensors + + //create the ladder support envelope + Box LadderSupportEnvelopeSolid(support_height/2.0, support_width/2.0, support_length/2.0); + Volume LadderSupportEnvelopeLogical(name + _toString( layer_id,"_SupEnvLogical_%02d"), LadderSupportEnvelopeSolid, air); + vxd.setVisAttributes(theDetector, "seeThrough", LadderSupportEnvelopeLogical); + + //create ladder support volume + Box LadderSupportSolid(support_thickness / 2.0 , support_width / 2.0 , support_length / 2.0); + Volume LadderSupportLogical(name + _toString( layer_id,"_SupLogical_%02d"), LadderSupportSolid, support_mat); + LadderSupportLogical.setVisAttributes(theDetector.visAttributes(supportVis)); + + //vxd.setVisAttributes(theDetector, sensVis, SensorLogical); + // vxd.setVisAttributes(theDetector, sensEnvVis, SensorBottomEnvelopeLogical); + // vxd.setVisAttributes(theDetector, ladderVis, LadderLogical); + + pv = LadderSupportEnvelopeLogical.placeVolume(LadderSupportLogical); + pv = LadderLogical.placeVolume(LadderSupportEnvelopeLogical); + + for(int i = 0; i < n_ladders; i++){ + std::stringstream ladder_enum; + ladder_enum << "vxt_ladder_" << layer_id << "_" << i; + DetElement ladderDE(layerDE, ladder_enum.str(), x_det.id()); + std::cout << "start building " << ladder_enum.str() << ":" << endl; + + //====== create the meassurement surface =================== + dd4hep::rec::Vector3D o(0,0,0); + dd4hep::rec::Vector3D u( 0., 0., 1.); + dd4hep::rec::Vector3D v( 0., 1., 0.); + dd4hep::rec::Vector3D n( 1., 0., 0.); + double inner_thick_top = sensor_thickness/2.0; + double outer_thick_top = support_height/2.0 + flex_thickness + sensor_thickness/2.0; + double inner_thick_bottom = support_height/2.0 + flex_thickness + sensor_thickness/2.0; + double outer_thick_bottom = sensor_thickness/2.0; + dd4hep::rec::VolPlane surfTop( SensorLogical , + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_top, outer_thick_top , u,v,n,o ) ; + dd4hep::rec::VolPlane surfBottom( SensorLogical , + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_bottom, outer_thick_bottom, u,v,n,o ) ; + + + for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ + std::stringstream topsensor_str; + std::stringstream bottomsensor_str; + topsensor_str << ladder_enum.str() << "_top_" << isensor; + // std::cout << "\tstart building " << topsensor_str.str() << ":" << endl; + bottomsensor_str << ladder_enum.str() << "_bottom_" << isensor; + // std::cout << "\tstart building " << bottomsensor_str.str() << ":" << endl; + DetElement topsensorDE(ladderDE, topsensor_str.str(), x_det.id()); + DetElement bottomsensorDE(ladderDE, bottomsensor_str.str(), x_det.id()); + topsensorDE.setPlacement(TopSensor_pv[isensor]); + volSurfaceList(topsensorDE)->push_back(surfTop); + // std::cout << "\t" << topsensor_str.str() << " done." << endl; + bottomsensorDE.setPlacement(BottomSensor_pv[isensor]); + // std::cout << "\t" << bottomsensor_str.str() << " done." << endl; + volSurfaceList(bottomsensorDE)->push_back(surfBottom); + } + Transform3D tr (RotationZYX(ladder_dphi*i,0.,0.),Position(ladder_radius*cos(ladder_phi0+ladder_dphi*i), ladder_radius*sin(ladder_phi0+ladder_dphi*i), 0.)); + pv = layer_assembly.placeVolume(LadderLogical,tr); + pv.addPhysVolID("module", i ) ; + ladderDE.setPlacement(pv); + std::cout << ladder_enum.str() << " done." << endl; + if(i==0) std::cout << "xy=" << ladder_radius*cos(ladder_phi0) << " " << ladder_radius*sin(ladder_phi0) << std::endl; + } + + // package the reconstruction data + dd4hep::rec::ZPlanarData::LayerLayout topLayer; + dd4hep::rec::ZPlanarData::LayerLayout bottomLayer; + + topLayer.ladderNumber = n_ladders; + topLayer.phi0 = 0.; + topLayer.sensorsPerLadder = n_sensors_per_side; + topLayer.lengthSensor = sensor_active_len; + topLayer.distanceSupport = sensitive_radius; + topLayer.thicknessSupport = support_thickness / 2.0; + topLayer.offsetSupport = -ladder_offset; + topLayer.widthSupport = support_width; + topLayer.zHalfSupport = support_length / 2.0; + topLayer.distanceSensitive = sensitive_radius + support_height / 2.0 + flex_thickness; + topLayer.thicknessSensitive = sensor_thickness; + topLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); + topLayer.widthSensitive = sensor_active_width; + topLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0; + + bottomLayer.ladderNumber = n_ladders; + bottomLayer.phi0 = 0.; + bottomLayer.sensorsPerLadder = n_sensors_per_side; + bottomLayer.lengthSensor = sensor_active_len; + bottomLayer.distanceSupport = sensitive_radius - support_height / 2.0 - flex_thickness; + bottomLayer.thicknessSupport = support_thickness / 2.0; + bottomLayer.offsetSupport = -ladder_offset; + bottomLayer.widthSupport = support_width; + bottomLayer.zHalfSupport = support_length / 2.0; + bottomLayer.distanceSensitive = sensitive_radius - support_height / 2.0 - sensor_thickness - flex_thickness; + bottomLayer.thicknessSensitive = sensor_thickness; + bottomLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); + bottomLayer.widthSensitive = sensor_active_width; + bottomLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0; + + zPlanarData->layers.push_back(bottomLayer); + zPlanarData->layers.push_back(topLayer); + } + std::cout << (*zPlanarData) << endl; + vxd.addExtension< ZPlanarData >(zPlanarData); + if ( x_det.hasAttr(_U(combineHits)) ) { + vxd.setCombineHits(x_det.attr<bool>(_U(combineHits)),sens); + } + std::cout << "vxd done." << endl; + return vxd; +} +DECLARE_DETELEMENT(SiTrackerStaggeredLadder_v01,create_element) diff --git a/Reconstruction/RecGenfitAlg/CMakeLists.txt b/Reconstruction/RecGenfitAlg/CMakeLists.txt index 6136171a0361d29af7c4541b92b1702bebb1b399..f7c5b6961e7e17e6bac3bfaabbe43e6efe07f8bf 100644 --- a/Reconstruction/RecGenfitAlg/CMakeLists.txt +++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt @@ -27,6 +27,7 @@ gaudi_add_module(RecGenfitAlg ) target_include_directories(RecGenfitAlg PUBLIC + ${LCIO_INCLUDE_DIRS} $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>/include $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}> ) diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp index 7e4dda5ce5b68c3bd6f9b03605b6b980e98b1654..df95acc6cbe9cfa3fe53902b659792f367b5c558 100644 --- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp @@ -32,6 +32,8 @@ //#include "Tools/KiTrackMarlinCEDTools.h" #include "Tools/FTDHelixFitter.h" +#include <TStopwatch.h> + using namespace MarlinTrk ; // Used to fedine the quality of the track output collection @@ -64,6 +66,23 @@ StatusCode ForwardTrackingAlg::initialize(){ // can be printed. As this is mainly used for debugging it is not a steerable parameter. //if( _useCED )MarlinCED::init(this) ; //CED + if(m_dumpTime){ + NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); + if ( !nt1 ) { + m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); + if ( 0 != m_tuple ) { + m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); + } + else { + fatal() << "Cannot book MyTuples/Time"+name() <<endmsg; + return StatusCode::FAILURE; + } + } + else{ + m_tuple = nt1; + } + } + // Now set min and max values for all the criteria for( unsigned i=0; i < _criteriaNames.size(); i++ ){ std::vector< float > emptyVec; @@ -181,6 +200,8 @@ StatusCode ForwardTrackingAlg::initialize(){ StatusCode ForwardTrackingAlg::execute(){ debug() << " processing event number " << _nEvt << endmsg; + auto stopwatch = TStopwatch(); + auto trkCol = _outColHdl.createAndPut(); ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // // @@ -683,6 +704,12 @@ StatusCode ForwardTrackingAlg::execute(){ //if( _useCED ) MarlinCED::draw(this); _nEvt ++ ; + + if(m_dumpTime&&m_tuple){ + m_timeTotal = stopwatch.RealTime()*1000; + m_tuple->write(); + } + return StatusCode::SUCCESS; } diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h index 56de6af9f2bb3ea94b24554b3139a10b658b0600..bd7a5173be02214e734579e3d4c025f8006ec1a2 100644 --- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h +++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h @@ -21,6 +21,8 @@ #include "Criteria/Criteria.h" #include "ILDImpl/SectorSystemFTD.h" +#include "GaudiKernel/NTuple.h" + using namespace KiTrack; using namespace KiTrackMarlin; @@ -234,7 +236,10 @@ class ForwardTrackingAlg : public GaudiAlgorithm { int _nRun ; int _nEvt ; - + + NTuple::Tuple* m_tuple; + NTuple::Item<float> m_timeTotal; + /** B field in z direction */ double _Bz; @@ -255,6 +260,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm { Gaudi::Property<std::vector<std::string> > _criteriaNames{this, "Criteria", Criteria::getAllCriteriaNamesVec()}; Gaudi::Property<std::vector<float> > _critMinimaInit{this, "CriteriaMin", {} }; Gaudi::Property<std::vector<float> > _critMaximaInit{this, "CriteriaMax", {} }; + Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true}; std::map<std::string, std::vector<float> > _critMinima; std::map<std::string, std::vector<float> > _critMaxima; diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp index 25589fad93f0265746c4fddcfb3fe71055b7fa24..f9b6a5466aba7cf08527d77d2f3b6c5d429ba3d8 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp @@ -34,6 +34,7 @@ #include "TrackSystemSvc/HelixFit.h" #include "TrackSystemSvc/IMarlinTrack.h" +#include <TStopwatch.h> //#include "TrackSystemSvc/MarlinTrkDiagnostics.h" //#ifdef MARLINTRK_DIAGNOSTICS_ON //#include "TrackSystemSvc/DiagnosticsController.h" @@ -104,7 +105,23 @@ StatusCode SiliconTrackingAlg::initialize() { _nRun = -1 ; _nEvt = 0 ; //printParameters() ; - + if(m_dumpTime){ + NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); + if ( !nt1 ) { + m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); + if ( 0 != m_tuple ) { + m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); + m_tuple->addItem ("timeKalman", m_timeKalman ).ignore(); + } + else { + fatal() << "Cannot book MyTuples/Time"+name() <<endmsg; + return StatusCode::FAILURE; + } + } + else{ + m_tuple = nt1; + } + } // set up the geometery needed by KalTest //FIXME: for now do KalTest only - make this a steering parameter to use other fitters auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc"); @@ -157,6 +174,7 @@ StatusCode SiliconTrackingAlg::initialize() { } StatusCode SiliconTrackingAlg::execute(){ + auto stopwatch = TStopwatch(); Navigation::Instance()->Initialize(); //_current_event = evt; //_allHits.reserve(1000); @@ -355,6 +373,10 @@ StatusCode SiliconTrackingAlg::execute(){ CleanUp(); debug() << "Event is done " << endmsg; _nEvt++; + if(m_dumpTime&&m_tuple){ + m_timeTotal = stopwatch.RealTime()*1000; + m_tuple->write(); + } return StatusCode::SUCCESS; } @@ -907,7 +929,7 @@ StatusCode SiliconTrackingAlg::finalize(){ //delete _trksystem ; _trksystem = 0; //delete _histos ; _histos = 0; info() << "Processed " << _nEvt << " events " << endmsg; - info() << lcio::ILDCellID0::encoder_string << " " << UTIL::ILDCellID0::encoder_string << endmsg; + return GaudiAlgorithm::finalize(); } @@ -2545,6 +2567,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { float pyTot = 0.; float pzTot = 0.; debug() << "Total " << nTracks << " candidate tracks will be dealed" << endmsg; + auto stopwatch = TStopwatch(); for (int iTrk=0;iTrk<nTracks;++iTrk) { TrackExtended * trackAR = _trackImplVec[iTrk]; @@ -2698,7 +2721,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { // setup initial dummy covariance matrix //std::vector<float> covMatrix; //covMatrix.resize(15); - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; @@ -2886,7 +2909,8 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { } } } - + if(m_dumpTime&&m_tuple) m_timeKalman = stopwatch.RealTime()*1000; + debug() << " -> run " << _nRun << " event " << _nEvt << endmsg; debug() << "Number of Si tracks = " << nSiSegments << endmsg; debug() << "Total 4-momentum of Si Tracks : E = " << std::setprecision(7) << eTot diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h index 7f584b264d3fe274bd7a793f3ac6ea0081249e28..e8e915b56530c595f1ed01040c2d5d6421af4bfb 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h @@ -27,6 +27,7 @@ #include <UTIL/BitField64.h> #include <UTIL/ILDConf.h> +#include "GaudiKernel/NTuple.h" //using namespace edm4hep ; namespace gear{ @@ -252,7 +253,7 @@ class SiliconTrackingAlg : public GaudiAlgorithm { Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true}; Gaudi::Property<bool> _SmoothOn{this, "SmoothOn", true}; Gaudi::Property<float> _helix_max_r{this, "HelixMaxR", 2000.}; - + Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true}; //std::vector<int> _colours; /** helper function to get collection using try catch block */ @@ -300,7 +301,11 @@ class SiliconTrackingAlg : public GaudiAlgorithm { std::vector<TrackerHitExtendedVec> _sectors; std::vector<TrackerHitExtendedVec> _sectorsFTD; - + + NTuple::Tuple* m_tuple; + NTuple::Item<float> m_timeTotal; + NTuple::Item<float> m_timeKalman; + /** * A helper class to allow good code readability by accessing tracks with N hits. * As the smalest valid track contains three hits, but the first index in a vector is 0, diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp index d918cf3831636acb238ac87a4ea0e7064e64f910..d4741f2eb165767d51fdaf7afee0612dacf942b7 100644 --- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp +++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp @@ -15,6 +15,8 @@ #include "TrackSystemSvc/MarlinTrkUtils.h" +#include <TStopwatch.h> + using namespace KiTrack; DECLARE_COMPONENT(TrackSubsetAlg) @@ -39,6 +41,23 @@ StatusCode TrackSubsetAlg::initialize() { _nRun = 0 ; _nEvt = 0 ; + if(m_dumpTime){ + NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); + if ( !nt1 ) { + m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); + if ( 0 != m_tuple ) { + m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); + } + else { + fatal() << "Cannot book MyTuples/Time"+name() <<endmsg; + return StatusCode::FAILURE; + } + } + else{ + m_tuple = nt1; + } + } + for(unsigned i=0; i<_trackInputColNames.size(); i++){ _inTrackColHdls.push_back(new DataHandle<edm4hep::TrackCollection> (_trackInputColNames[i], Gaudi::DataHandle::Reader, this)); } @@ -96,6 +115,8 @@ StatusCode TrackSubsetAlg::finalize(){ } StatusCode TrackSubsetAlg::execute(){ + auto stopwatch = TStopwatch(); + std::vector<edm4hep::Track> tracks; auto trkCol = _outColHdl.createAndPut(); @@ -214,7 +235,7 @@ StatusCode TrackSubsetAlg::execute(){ trackerHits.push_back(Navigation::Instance()->GetTrackerHit(trackerHitsObj[i].getObjectID())); } // setup initial dummy covariance matrix - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; } @@ -331,6 +352,12 @@ StatusCode TrackSubsetAlg::execute(){ Navigation::Instance()->Initialize(); _nEvt ++ ; + + if(m_dumpTime&&m_tuple){ + m_timeTotal = stopwatch.RealTime()*1000; + m_tuple->write(); + } + return StatusCode::SUCCESS; } diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h index 952d08fa2607dee095444f3dcf0563901c28dc63..1023eba976e861461174d93a7158294328ef4854 100644 --- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h +++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h @@ -11,6 +11,7 @@ #include "Math/ProbFunc.h" +#include "GaudiKernel/NTuple.h" /** Processor that takes tracks from multiple sources and outputs them (or modified versions, or a subset of them) * as one track collection. * @@ -74,11 +75,15 @@ class TrackSubsetAlg : public GaudiAlgorithm { Gaudi::Property<float> _initialTrackError_tanL{this, "InitialTrackErrorTanL",1e2}; Gaudi::Property<double> _maxChi2PerHit{this, "MaxChi2PerHit", 1e2}; Gaudi::Property<double> _omega{this, "Omega", 0.75}; + Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true}; float _bField; int _nRun ; int _nEvt ; + + NTuple::Tuple* m_tuple; + NTuple::Item<float> m_timeTotal; }; /** A functor to return whether two tracks are compatible: The criterion is if the share a TrackerHit or more */ diff --git a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h index c12f55174033068d259fa520eed26bd16fd47eef..6235b3fb25e6a06ed8d6d8ceee9e8127ad370f71 100644 --- a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h +++ b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h @@ -27,9 +27,9 @@ inline edm4hep::TrackState getTrackStateAt(edm4hep::Track track, int location) { return edm4hep::TrackState(); } -inline std::array<float,15> getCovMatrix(const edm4hep::Track &track) { +inline decltype(edm4hep::TrackState::covMatrix) getCovMatrix(const edm4hep::Track &track) { if(track.trackStates_size()>0) return track.getTrackStates(0).covMatrix; - std::array<float,15> dummy{}; + decltype(edm4hep::TrackState::covMatrix) dummy{}; return dummy; } inline float getTanLambda(const edm4hep::Track &track) { diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp index 9d57f09cae8a4f8a25bd1a811b571638ec15cb30..9ecee56418ae4387e73572d73b42b757a92bd2b2 100644 --- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp +++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp @@ -1238,31 +1238,35 @@ start: bool reverse_order = ( std::abs( hf->first->pos.z() ) > std::abs( hb->first->pos.z()) + 3. ) ; unsigned nHit = 0 ; - + int code = 0; if( reverse_order ){ //std::cout << "It is true order" << std::endl; for( CluTrack::reverse_iterator it=clu->rbegin() ; it != clu->rend() ; ++it){ edm4hep::TrackerHit ph = (*it)->first->edm4hepHit; trk->addHit(ph) ; ++nHit ; - //std::cout << " hit added " << (*it)->first->edm4hepHit << std::endl ; + //std::cout << " hit added " << (*it)->first->edm4hepHit.id() << std::endl ; } - trk->initialise( MarlinTrk::IMarlinTrack::forward ) ; + code = trk->initialise( MarlinTrk::IMarlinTrack::forward ) ; } else { //std::cout << "It is reverse order" << std::endl; for( CluTrack::iterator it=clu->begin() ; it != clu->end() ; ++it){ edm4hep::TrackerHit ph = (*it)->first->edm4hepHit; - trk->addHit(ph) ; + if( trk->addHit(ph) == MarlinTrk::IMarlinTrack::success ){ + //std::cout << " hit added " << (*it)->first->edm4hepHit.id() << std::endl; + } + else{ + //std::cout << " hit not added " << (*it)->first->edm4hepHit.id() << std::endl; + } ++nHit ; - //std::cout << " hit added "<< (*it)->first->edm4hepHit << std::endl ; } - trk->initialise( MarlinTrk::IMarlinTrack::backward ) ; + code = trk->initialise( MarlinTrk::IMarlinTrack::backward ) ; } - int code = trk->fit( maxChi2 ) ; + if( code != MarlinTrk::IMarlinTrack::error ) code = trk->fit( maxChi2 ) ; // for (auto hit : hitsInFit) std::cout << hit.first << std::endl; if( code != MarlinTrk::IMarlinTrack::success ){ @@ -1270,7 +1274,8 @@ start: std::cout << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IMarlinTrkFitter : problem fitting track " << " error code : " << MarlinTrk::errorCode( code ) << std::endl ; - + delete trk; + return 0; } @@ -1361,7 +1366,7 @@ start: tsBase.Z0 = 0; tsBase.tanLambda = 0; tsBase.referencePoint = edm4hep::Vector3f(0,0,0); - tsBase.covMatrix = std::array<float, 15>{}; + tsBase.covMatrix = decltype(edm4hep::TrackState::covMatrix){}; edm4hep::TrackState tsIP(tsBase); edm4hep::TrackState tsFH(tsBase); edm4hep::TrackState tsLH(tsBase); diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h index 46247d0b48c36b6975da0ba3411f27f9363e99d8..6ace630a7dcf6bdc552426f519f62c51a992da1a 100644 --- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h +++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h @@ -467,11 +467,12 @@ namespace clupatra_new{ // FIXME Mingrui debug // streamlog_out( DEBUG3 ) << " -- extrapolate TrackState : " << lcshort( ts ) << std::endl ; - - edm4hep::TrackerHit ht = trk.getTrackerHits(0); - //need to add a dummy hit to the track - mTrk->addHit( ht ) ; // is this the right hit ?????????? - + // fucd: getTrackerHits(0) is possible to miss ILDVTrackHit + for(int ih=0;ih<nHit;ih++){ + edm4hep::TrackerHit ht = trk.getTrackerHits(ih); + //need to add a dummy hit to the track + if(mTrk->addHit( ht ) == MarlinTrk::IMarlinTrack::success) break; // is this the right hit ?????????? + } mTrk->initialise( ts , _b , MarlinTrk::IMarlinTrack::backward ) ; double deltaChi ; diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp index f5a986328dd86bf39acd24f0bfad6607c11d3fec..08e1fb022fdfc240dd127957085fb954ca2edf54 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp @@ -53,6 +53,8 @@ #include <vector> #include <bitset> +#include <TStopwatch.h> + typedef std::vector<edm4hep::TrackerHit> TrackerHitVec; using namespace edm4hep ; @@ -143,7 +145,24 @@ StatusCode FullLDCTrackingAlg::initialize() { PI = acos(-1.); PIOVER2 = 0.5*PI; TWOPI = 2*PI; - + + if(m_dumpTime){ + NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); + if ( !nt1 ) { + m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); + if ( 0 != m_tuple ) { + m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); + m_tuple->addItem ("timeKalman", m_timeKalman ).ignore(); + } + else { + fatal() << "Cannot book MyTuples/Time"+name() <<endmsg; + return StatusCode::FAILURE; + } + } + else{ + m_tuple = nt1; + } + } // set up the geometery needed by KalTest //FIXME: for now do KalTest only - make this a steering parameter to use other fitters auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc"); @@ -181,6 +200,8 @@ StatusCode FullLDCTrackingAlg::execute() { // debug() << endmsg; info() << "FullLDCTrackingAlg -> run = " << 0/*evt->getRunNumber()*/ << " event = " << _nEvt << endmsg; + + auto stopwatch = TStopwatch(); // debug() << endmsg; auto outCol = _OutputTrackColHdl.createAndPut(); @@ -222,6 +243,11 @@ StatusCode FullLDCTrackingAlg::execute() { debug() << "Cleanup is done." << endmsg; _nEvt++; // getchar(); + if(m_tuple){ + m_timeTotal = stopwatch.RealTime()*1000; + m_tuple->write(); + } + debug() << "FullLDCTrackingAlg execute() finished" << endmsg; return StatusCode::SUCCESS; @@ -248,6 +274,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr //SJA:FIXME: So here we are going to do one final refit. This can certainly be optimised, but rather than worry about the mememory management right now lets make it work, and optimise it later ... debug() << "Total " << nTrkCand << " trkCand to deal" << endmsg; + auto stopwatch = TStopwatch(); for (int iTRK=0;iTRK<nTrkCand;++iTRK) { TrackExtended * trkCand = trkVec[iTRK]; @@ -288,7 +315,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr edm4hep::MutableTrack track;// = new edm4hep::Track; // setup initial dummy covariance matrix - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; @@ -385,9 +412,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr int error_code = 0; try { - error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit); - } catch (...) { // delete track; @@ -554,7 +579,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr debug() << " Add Track to final Collection: ID = " << track.id() << " for trkCand "<< trkCand << endmsg; } } - + if(m_tuple) m_timeKalman = stopwatch.RealTime()*1000; // streamlog_out(DEBUG5) << endmsg; debug() << "Number of accepted " << _OutputTrackColHdl.fullKey() << " = " << nTotTracks << endmsg; debug() << "Total 4-momentum of " << _OutputTrackColHdl.fullKey() << " : E = " << eTot @@ -1115,7 +1140,7 @@ void FullLDCTrackingAlg::prepareVectors() { //param[3] = getD0(tpcTrack); //param[4] = getZ0(tpcTrack); - std::array<float, 15> Cov = getCovMatrix(tpcTrack); + auto Cov = getCovMatrix(tpcTrack); int NC = int(Cov.size()); for (int ic=0;ic<NC;ic++) { cov[ic] = Cov[ic]; @@ -1179,7 +1204,7 @@ void FullLDCTrackingAlg::prepareVectors() { //param[3] = getD0(siTrack); //param[4] = getZ0(siTrack); - std::array<float, 15> Cov = getCovMatrix(siTrack); + auto Cov = getCovMatrix(siTrack); int NC = int(Cov.size()); for (int ic=0;ic<NC;ic++) { cov[ic] = Cov[ic]; @@ -1571,7 +1596,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac } // setup initial dummy covariance matrix - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; @@ -3411,7 +3436,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() { nonAssignedTPCHits.push_back(trkHitExt); } } - debug() << "AddNotAssignedHits : Number of Non Assigned TPC hits = " << nonAssignedTPCHits.size() << endmsg; + debug() << "AddNotAssignedHits : Number of Non Assigned TPC hits = " << nonAssignedTPCHits.size() << " distance cut = " << _distCutForTPCHits << endmsg; AssignTPCHitsToTracks(nonAssignedTPCHits, _distCutForTPCHits); } @@ -3628,7 +3653,7 @@ void FullLDCTrackingAlg::AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, f pre_fit.location = 1/*lcio::TrackState::AtIP*/; // setup initial dummy covariance matrix - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; @@ -3854,6 +3879,7 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec, debug() << "AssignTPCHitsToTracks: Starting loop " << nTrk << " tracks and " << nHits << " hits" << endmsg; for (int iT=0;iT<nTrk;++iT) { // loop over all tracks + debug() << " track " << iT << endmsg; TrackExtended * foundTrack = _trkImplVec[iT]; int tanlambdaSign = std::signbit(foundTrack->getTanLambda());//we only care about positive or negative GroupTracks * group = foundTrack->getGroupTracks(); @@ -3872,7 +3898,7 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec, HelixClass helix; helix.Initialize_Canonical(phi0,d0,z0,omega,tanLambda,_bField); float OnePFivehalfPeriodZ = 1.5*fabs(acos(-1.)*tanLambda/omega); - + debug() << " OnePFivehalfPeriodZ = " << OnePFivehalfPeriodZ << endmsg; for (int iH=0;iH<nHits;++iH) { // loop over leftover TPC hits //check if the hit and the track or on the same side @@ -3884,9 +3910,10 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec, bool consider = DeltaStart <= OnePFivehalfPeriodZ; consider = consider || (DeltaEnd <= OnePFivehalfPeriodZ); consider = consider || ( (HitPositions[iH][2]>=startPointZ) && (HitPositions[iH][2]<=endPointZ) ); - + debug() << " hit " << iH << " " << consider << " DeltaStart = " << DeltaStart << " DeltaEnd = " << DeltaEnd << endmsg; if(consider){ float distance = helix.getDistanceToPoint(HitPositions[iH], minDistances[iH]); + debug() << " distance = " << distance << " cut = " << minDistances[iH] << endmsg; if (distance < minDistances[iH]) { minDistances[iH] = distance; tracksToAttach[iH] = foundTrack; @@ -3901,7 +3928,9 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec, if (tracksToAttach[iH]!=NULL) { tracksToAttach[iH]->addTrackerHitExtended(trkHitExt); trkHitExt->setTrackExtended( tracksToAttach[iH] ); - trkHitExt->setUsedInFit( false ); + //by fucd + //trkHitExt->setUsedInFit( false ); + trkHitExt->setUsedInFit( true ); } } @@ -4134,7 +4163,7 @@ void FullLDCTrackingAlg::AssignSiHitsToTracks(TrackerHitExtendedVec hitVec, pre_fit.location = 1/*lcio::TrackState::AtIP*/; // setup initial dummy covariance matrix - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h index 6822d60f4dd05007b757256eb2e57efe13435ecc..cfc804879d44a7e453d069e79abe261640494120 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h @@ -27,6 +27,7 @@ #include <UTIL/BitField64.h> #include <UTIL/ILDConf.h> +#include "GaudiKernel/NTuple.h" //using namespace edm4hep ; //class gear::GearMgr ; @@ -394,7 +395,7 @@ protected: Gaudi::Property<float> _vetoMergeMomentumCut{this, "VetoMergeMomentumCut", 2.5}; Gaudi::Property<float> _maxAllowedPercentageOfOutliersForTrackCombination{this, "MaxAllowedPercentageOfOutliersForTrackCombination", 0.3}; Gaudi::Property<int> _maxAllowedSiHitRejectionsForTrackCombination{this, "MaxAllowedSiHitRejectionsForTrackCombination", 2}; - + Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true}; //float _dPCutForForcedMerging; bool _runMarlinTrkDiagnostics = false; @@ -514,6 +515,10 @@ protected: DataHandle<edm4hep::TrackCollection> _SiTrackColHdl{"SiTracks", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackCollection> _OutputTrackColHdl{"MarlinTrkTracks", Gaudi::DataHandle::Writer, this}; + + NTuple::Tuple* m_tuple; + NTuple::Item<float> m_timeTotal; + NTuple::Item<float> m_timeKalman; }; #endif diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp index c87ebc39a212dfaa2ac67d6780f580aa6a17728e..afe625ae476eaf3afdcf32ac67149fc5276e193d 100644 --- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp +++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp @@ -458,8 +458,10 @@ void TruthTrackerAlg::getTrackStateFromMcParticle( trackState.Z0=helix.getZ0(); trackState.tanLambda=helix.getTanLambda(); trackState.referencePoint=helix.getReferencePoint(); - std::array<float,15> covMatrix; - for(int i=0;i<15;i++){covMatrix[i]=1.;}//FIXME + + decltype(trackState.covMatrix) covMatrix; + for(int i=0;i<covMatrix.size();i++){covMatrix[i]=999.;}//FIXME + trackState.covMatrix=covMatrix; getCircleFromPosMom(pos,mom,B[2]/dd4hep::tesla,mcParticle.getCharge(),m_helixRadius,m_helixXC,m_helixYC); diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp index 81c3611f63eafe0e0eb537b127b4cb44c0045a74..e96ae681694de49f30b430dedf1ec1ac9568f0e3 100644 --- a/Service/GearSvc/src/GearSvc.cpp +++ b/Service/GearSvc/src/GearSvc.cpp @@ -12,11 +12,13 @@ #include "gearimpl/FixedPadSizeDiskLayout.h" #include "gearimpl/CalorimeterParametersImpl.h" #include "gearimpl/SimpleMaterialImpl.h" +#include "gear/VXDLayerLayout.h" #include "gearxml/tinyxml.h" #include "DD4hep/Detector.h" #include "DD4hep/DetElement.h" #include "DDRec/DetectorData.h" +#include "DDRec/MaterialManager.h" #include "DD4hep/DD4hepUnits.h" #include "CLHEP/Units/SystemOfUnits.h" @@ -194,6 +196,49 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){ extensionDataValid = false; info() << e.what() << " " << vxdData << endmsg; } + if(vxdData){ + int vxdType = gear::ZPlanarParameters::CMOS; + gear::ZPlanarParametersImpl* gearVXD = new gear::ZPlanarParametersImpl( vxdType, vxdData->rInnerShell/dd4hep::mm, vxdData->rOuterShell/dd4hep::mm, + vxdData->zHalfShell/dd4hep::mm , vxdData->gapShell/dd4hep::mm , 0. ) ; + for(unsigned i=0,n=vxdData->layers.size() ; i<n; ++i){ + const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[i]; + // FIXME set rad lengths to 0 -> need to get from dd4hep .... + gearVXD->addLayer(l.ladderNumber, l.phi0, + l.distanceSupport/dd4hep::mm, l.offsetSupport/dd4hep::mm, l.thicknessSupport/dd4hep::mm, l.zHalfSupport/dd4hep::mm, l.widthSupport/dd4hep::mm, 0., + l.distanceSensitive/dd4hep::mm, l.offsetSensitive/dd4hep::mm, l.thicknessSensitive/dd4hep::mm, l.zHalfSensitive/dd4hep::mm, l.widthSensitive/dd4hep::mm, 0.); + } + m_gearMgr->setVXDParameters(gearVXD); + + dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ; + const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[0] ; + dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ; + dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ; + + std::cout << " ####### found materials between points : " << a << " and " << b << " : " ; + for( unsigned i=0,n=materials.size();i<n;++i){ + std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ; + } + std::cout << std::endl ; + std::cout << " averaged material : " << mat << std::endl ; + gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", mat.A(), mat.Z(), + mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), + mat.radiationLength()/dd4hep::mm, + mat.interactionLength()/dd4hep::mm); + m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); + + info() << vxdData->rInnerShell << " " << vxdData->rOuterShell << " " << vxdData->zHalfShell << " " << vxdData->gapShell << endmsg; + for(int i=0,n=vxdData->layers.size(); i<n; i++){ + const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = vxdData->layers[i]; + info() << i << ": " << thisLayer.ladderNumber << "," << thisLayer.phi0 << "," + << thisLayer.distanceSupport/dd4hep::mm << "," << thisLayer.offsetSupport/dd4hep::mm << "," << thisLayer.thicknessSupport/dd4hep::mm << "," + << thisLayer.zHalfSupport/dd4hep::mm << "," << thisLayer.widthSupport/dd4hep::mm << "," << "NULL," + << thisLayer.distanceSensitive/dd4hep::mm << "," << thisLayer.offsetSensitive/dd4hep::mm << "," << thisLayer.thicknessSensitive/dd4hep::mm << "," + << thisLayer.zHalfSensitive/dd4hep::mm << "," << thisLayer.widthSensitive/dd4hep::mm << ",NULL" << endmsg; + } + return sc; + } std::vector<helpLayer> helpSensitives; std::vector<helpLayer> helpLadders; @@ -575,30 +620,34 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){ gearParameters->setDoubleVal("CryostatAlHalfZ", dzSty+drSty); m_gearMgr->setGearParameters("VXDInfra", gearParameters); } + + dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ; + const gear::VXDLayerLayout& l = m_gearMgr->getVXDParameters().getVXDLayerLayout(); + dd4hep::rec::Vector3D a( l.getSensitiveDistance(0)+l.getSensitiveThickness(0), l.getPhi0(0), 0. , dd4hep::rec::Vector3D::cylindrical ) ; + dd4hep::rec::Vector3D b( l.getLadderDistance(0)+l.getLadderThickness(0), l.getPhi0(0), 0. , dd4hep::rec::Vector3D::cylindrical ) ; + const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ; + dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ; + + std::cout << " ####### found materials between points : " << a << " and " << b << " : " ; + for( unsigned i=0,n=materials.size();i<n;++i){ + std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ; + } + std::cout << std::endl ; + std::cout << " averaged material : " << mat << std::endl ; + gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", mat.A(), mat.Z(), + mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), + mat.radiationLength()/dd4hep::mm, + mat.interactionLength()/dd4hep::mm); + //effective A different with what in Mokka, fix them as Mokka's gear::SimpleMaterialImpl* VXDFoamShellMaterial_old = new gear::SimpleMaterialImpl("VXDFoamShellMaterial_old", 1.043890843e+01, 5.612886646e+00, 2.500000000e+01, 1.751650267e+04, 0); m_gearMgr->registerSimpleMaterial(VXDFoamShellMaterial_old); gear::SimpleMaterialImpl* VXDFoamShellMaterial = new gear::SimpleMaterialImpl("VXDFoamShellMaterial", styAeff, styZeff, styDensity*(CLHEP::g/CLHEP::cm3)/(CLHEP::kg/CLHEP::m3), styRadLen, styIntLen); m_gearMgr->registerSimpleMaterial(VXDFoamShellMaterial); - gear::SimpleMaterialImpl* VXDSupportMaterial_old = new gear::SimpleMaterialImpl("VXDSupportMaterial_old", 2.075865162e+01, 1.039383117e+01, 2.765900000e+02, 1.014262421e+03, 3.341388059e+03); - m_gearMgr->registerSimpleMaterial(VXDSupportMaterial_old); - gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", VXDSupportAeff, VXDSupportZeff, VXDSupportDensity, VXDSupportRadLen, VXDSupportIntLen); m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); - info() << "=====================from ZPlanarData==============================" << endmsg; - if(vxdData){ - info() << vxdData->rInnerShell << " " << vxdData->rOuterShell << " " << vxdData->zHalfShell << " " << vxdData->gapShell << endmsg; - const std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& layers = vxdData->layers; - for(int i=0;i<layers.size();i++){ - const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = layers[i]; - info() << i << ": " << thisLayer.ladderNumber << "," << thisLayer.phi0 << "," << thisLayer.distanceSupport << "," << thisLayer.offsetSupport << "," - << thisLayer.thicknessSupport << "," << thisLayer.zHalfSupport << "," << thisLayer.widthSupport << "," << "NULL," - << thisLayer.distanceSensitive << "," << thisLayer.offsetSensitive << "," << thisLayer.thicknessSensitive << "," << thisLayer.zHalfSensitive << "," - << thisLayer.widthSensitive << ",NULL" << endmsg; - } - } - info() << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg; - //info() << m_materials["VXDSupportMaterial"] << endmsg; + info() << "cryostat parameters: " << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg; + return sc; } @@ -853,7 +902,7 @@ StatusCode GearSvc::convertDC(dd4hep::DetElement& dc){ } info() << (*dcData) << endmsg; } - + debug() << dcData->maxRow << ": " << dcData->rMinReadout*CLHEP::cm << " " << dcData->rMaxReadout*CLHEP::cm << endmsg; // regard as TPCParameters, TODO: drift chamber parameters gear::TPCParametersImpl *tpcParameters = new gear::TPCParametersImpl(); gear::PadRowLayout2D *padLayout = new gear::FixedPadSizeDiskLayout(dcData->rMinReadout*CLHEP::cm, dcData->rMaxReadout*CLHEP::cm, diff --git a/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h b/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h index 68498fe0dd6ff433882733ec36da023602edae25..70303ae34a3c9547e2cfbe359e7cc651aace6287 100644 --- a/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h +++ b/Service/TrackSystemSvc/include/TrackSystemSvc/MarlinTrkUtils.h @@ -55,7 +55,7 @@ namespace MarlinTrk{ std::vector<edm4hep::TrackerHit>& hit_list, edm4hep::MutableTrack* track, bool fit_backwards, - const std::array<float,15>& initial_cov_for_prefit, + const decltype(edm4hep::TrackState::covMatrix)& initial_cov_for_prefit, float bfield_z, double maxChi2Increment=DBL_MAX); diff --git a/Service/TrackSystemSvc/src/LCIOTrackPropagators.cc b/Service/TrackSystemSvc/src/LCIOTrackPropagators.cc index 22fb8665f77fe8dfb162997685f86224514d2218..a5ab0265dd3d21e32dcab7c48eb543da39d6968e 100644 --- a/Service/TrackSystemSvc/src/LCIOTrackPropagators.cc +++ b/Service/TrackSystemSvc/src/LCIOTrackPropagators.cc @@ -100,7 +100,7 @@ namespace LCIOTrackPropagators{ CLHEP::HepSymMatrix covPrime = cov0.similarity(propagatorMatrix); - std::array<float,15> cov; + decltype(edm4hep::TrackState::covMatrix) cov; icov = 0 ; diff --git a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc index 7979379267aa921069505248c3a65aa6ebaffd44..ffedd64edb8558338e69a028ddcdcbf8021d42e1 100644 --- a/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc +++ b/Service/TrackSystemSvc/src/MarlinKalTestTrack.cc @@ -376,7 +376,7 @@ namespace MarlinTrk { bfield_z ); TMatrixD cov(5,5) ; - std::array<float, 15> covLCIO = ts.covMatrix; + auto covLCIO = ts.covMatrix; cov( 0 , 0 ) = covLCIO[ 0] ; // d0, d0 cov( 0 , 1 ) = - covLCIO[ 1] ; // d0, phi @@ -1599,7 +1599,7 @@ namespace MarlinTrk { Double_t cpa = helix.GetKappa(); double alpha = omega / cpa ; // conversion factor for omega (1/R) to kappa (1/Pt) - std::array<float, 15> covLCIO; + decltype(ts.covMatrix) covLCIO; covLCIO[ 0] = covK( 0 , 0 ) ; // d0, d0 covLCIO[ 1] = - covK( 1 , 0 ) ; // phi0, d0 diff --git a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc index 3c2dead214e4fda49a0730fd13b8ee2843e19178..ffa0f7f9395accd360f1efce633308447873d258 100644 --- a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc +++ b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc @@ -94,7 +94,7 @@ namespace MarlinTrk { 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 std::array<float,15>& initial_cov_for_prefit, float bfield_z, double maxChi2Increment){ + 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){ /////////////////////////////////////////////////////// // check inputs diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp index ff20af31bc40cd2d24304e83e2c9bbadffc6a171..fb2c1126b79293482e3c66c6ebb3ffe9681302d7 100644 --- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp +++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp @@ -76,6 +76,14 @@ AnExampleDetElemTool::getLV() { void AnExampleDetElemTool::ConstructSDandField() { + + // DEBUG ONLY: turn off all the SD. + if (not m_SD_enabled) { + warning() << "All the Sensitive Detectors will be disabled by default. " << endmsg; + return; + } + + // // Construct SD using DD4hep. // Refer to FCCSW/Detector/DetComponents/src/ diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h index 0fd9c35627c5dfc56e4d09fa2370bfc79ce827e0..1b3a923f1f317a3c5e35b9d5d54eb30c84a0b097 100644 --- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h +++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h @@ -31,6 +31,8 @@ private: // DD4hep XML compact file path Gaudi::Property<std::string> m_dd4hep_xmls{this, "detxml"}; + Gaudi::Property<bool> m_SD_enabled{this, "SDenabled", true}; + SmartIF<IGeomSvc> m_geosvc; ToolHandle<ISensDetTool> m_calo_sdtool; ToolHandle<ISensDetTool> m_driftchamber_sdtool; diff --git a/Simulation/DetSimSD/include/DetSimSD/CaloSensitiveDetector.h b/Simulation/DetSimSD/include/DetSimSD/CaloSensitiveDetector.h index db161be4269e9ba94c071a50f647b5a99b38ec19..4f2b6182c431cf3249591ae06f435bec7240b38d 100644 --- a/Simulation/DetSimSD/include/DetSimSD/CaloSensitiveDetector.h +++ b/Simulation/DetSimSD/include/DetSimSD/CaloSensitiveDetector.h @@ -23,6 +23,7 @@ public: virtual void Initialize(G4HCofThisEvent* HCE); virtual G4bool ProcessHits(G4Step* step,G4TouchableHistory* history); virtual void EndOfEvent(G4HCofThisEvent* HCE); + void ApplyBirksLaw(){m_applyBirksLaw = true;}; protected: CalorimeterHit* find(const HitCollection*, const dd4hep::sim::HitCompare<CalorimeterHit>&); @@ -31,7 +32,8 @@ protected: HitCollection* m_hc; std::map<unsigned long, CalorimeterHit*> m_hitMap; - bool m_isMergeEnabled; + bool m_isMergeEnabled = false; + bool m_applyBirksLaw = false; }; diff --git a/Simulation/DetSimSD/src/CaloSensitiveDetector.cpp b/Simulation/DetSimSD/src/CaloSensitiveDetector.cpp index 11fb8cddabccf91713e068f0a4396772079ce2c1..950c4c45f91c1084f21fc95982fb533ae6b7a596 100644 --- a/Simulation/DetSimSD/src/CaloSensitiveDetector.cpp +++ b/Simulation/DetSimSD/src/CaloSensitiveDetector.cpp @@ -33,6 +33,7 @@ CaloSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { // std::cout << "CaloSensitiveDetector::ProcessHits" << std::endl; dd4hep::sim::Geant4StepHandler h(step); + if(m_applyBirksLaw) h.doApplyBirksLaw(); dd4hep::Position pos = 0.5 * (h.prePos() + h.postPos()); HitContribution contrib = dd4hep::sim::Geant4Hit::extractContribution(step); const std::string& name = GetName(); @@ -62,9 +63,9 @@ CaloSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { m_hc->insert(hit); } hit->truth.push_back(contrib); - hit->energyDeposit += contrib.deposit; - - + //hit->energyDeposit += contrib.deposit; + hit->energyDeposit += h.totalEnergy(); + //std::cout << "Apply Birk law: before = " << contrib.deposit << " after = " << h.totalEnergy() << std::endl; return true; } diff --git a/Simulation/DetSimSD/src/CalorimeterSensDetTool.cpp b/Simulation/DetSimSD/src/CalorimeterSensDetTool.cpp index 4219802f2f0df1bdc46be9552fd9661b3376295e..12a5f3f9586884efa9f5f3ee95f5aa9a82281362 100644 --- a/Simulation/DetSimSD/src/CalorimeterSensDetTool.cpp +++ b/Simulation/DetSimSD/src/CalorimeterSensDetTool.cpp @@ -42,8 +42,18 @@ CalorimeterSensDetTool::createSD(const std::string& name) { break; } } - G4VSensitiveDetector* sd = new CaloSensitiveDetector(name, *dd4hep_geo, is_merge_enabled); - debug() << name << " set to merge true/false = " << is_merge_enabled << endmsg; + + CaloSensitiveDetector* sd = new CaloSensitiveDetector(name, *dd4hep_geo, is_merge_enabled); + warning() << name << " set to merge true/false = " << is_merge_enabled << endmsg; + + for(auto cal_name : m_listCalsApplyBirks){ + if(cal_name==name){ + info() << name << " will apply Birks law" << endmsg; + sd->ApplyBirksLaw(); + break; + } + } + return sd; } diff --git a/Simulation/DetSimSD/src/CalorimeterSensDetTool.h b/Simulation/DetSimSD/src/CalorimeterSensDetTool.h index 549de18d10ee548d6ba2daa3f2734d65f3b156ee..69bab1d4e5a5bde38b087ce14d54c8db1969dd96 100644 --- a/Simulation/DetSimSD/src/CalorimeterSensDetTool.h +++ b/Simulation/DetSimSD/src/CalorimeterSensDetTool.h @@ -30,6 +30,7 @@ private: SmartIF<IGeomSvc> m_geosvc; Gaudi::Property<std::vector<std::string> > m_listCalsMergeDisable{this, "CalNamesMergeDisable", {}}; + Gaudi::Property<std::vector<std::string> > m_listCalsApplyBirks{this, "CalNamesApplyBirks", {}}; }; #endif diff --git a/Utilities/KiTrack/src/Tools/Fitter.cc b/Utilities/KiTrack/src/Tools/Fitter.cc index b7a0c156d760041719a7ddf0b095fa96a6265057..436f7ed222cecf80b75bec6cf5aa284751e11710 100644 --- a/Utilities/KiTrack/src/Tools/Fitter.cc +++ b/Utilities/KiTrack/src/Tools/Fitter.cc @@ -166,7 +166,7 @@ void Fitter::fitVXD(){ /**********************************************************************************************/ /* Create a TrackStateImpl from the helix values and use it to initalise the fit */ /**********************************************************************************************/ - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; @@ -185,6 +185,7 @@ void Fitter::fitVXD(){ helixTrack.getOmega(), helixTrack.getZ0(), helixTrack.getTanLambda(), + 0.f, // dummy value for time referencePoint, covMatrix}; @@ -349,7 +350,7 @@ void Fitter::fit(){ /**********************************************************************************************/ /* Create a TrackStateImpl from the helix values and use it to initalise the fit */ /**********************************************************************************************/ - std::array<float,15> covMatrix; + decltype(edm4hep::TrackState::covMatrix) covMatrix; for (unsigned icov = 0; icov<covMatrix.size(); ++icov) { covMatrix[icov] = 0; } @@ -365,7 +366,8 @@ void Fitter::fit(){ helixTrack.getPhi0(), helixTrack.getOmega(), helixTrack.getZ0(), - helixTrack.getTanLambda(), + helixTrack.getTanLambda(), + 0.f, // dummy value for time referencePoint, covMatrix}; diff --git a/build.sh b/build.sh index 8f8ad7e8c2a764f777be25a131acd812214f069c..337287fe081bdbf364c5fa2b856dc0a4534aa9df 100755 --- a/build.sh +++ b/build.sh @@ -92,7 +92,7 @@ function run-make() { # The current default platform lcg_platform=x86_64-centos7-gcc8-opt -lcg_version=101.0.0 +lcg_version=101.0.1 bldtool=${CEPCSW_BLDTOOL} # make, ninja # set in env var diff --git a/setup.sh b/setup.sh index 6b973779f7ef25ec5b49b4a97a0c9e4714217961..a75d5f39d2d3fca4d0501c970c74648b69229014 100644 --- a/setup.sh +++ b/setup.sh @@ -49,7 +49,7 @@ function setup-external() { # CEPCSW_LCG_VERSION=${1}; shift if [ -z "$CEPCSW_LCG_VERSION" ]; then - CEPCSW_LCG_VERSION=101.0.0 + CEPCSW_LCG_VERSION=101.0.1 fi export CEPCSW_LCG_VERSION