diff --git a/Detector/DetCRD/CMakeLists.txt b/Detector/DetCRD/CMakeLists.txt index 4a92fed45c70d5483398dc43e4e8abecd79ef302..d3f0b06954b5d3a547cb13495dfb52201ebf0818 100644 --- a/Detector/DetCRD/CMakeLists.txt +++ b/Detector/DetCRD/CMakeLists.txt @@ -17,6 +17,8 @@ find_package(ROOT COMPONENTS MathCore GenVector Geom REQUIRED) gaudi_add_module(DetCRD SOURCES src/Calorimeter/CRDEcal_v01.cpp + src/Calorimeter/RotatedPolyhedraBarrelCalorimeter_v01_geo.cpp + src/Calorimeter/RotatedCrystalCalorimeter_v01_geo.cpp src/Other/CRDBeamPipe_v01_geo.cpp src/Tracker/SiTrackerSkewRing_v01_geo.cpp diff --git a/Detector/DetCRD/compact/CRD_common_v01/Ecal_Rotated_Crystal_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/Ecal_Rotated_Crystal_v01_01.xml new file mode 100644 index 0000000000000000000000000000000000000000..eef2ba3cfdbe9e3c0dba6ad0c913ed870a3f41ef --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v01/Ecal_Rotated_Crystal_v01_01.xml @@ -0,0 +1,36 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd> + <define> + <constant name="Ecal_crystal_y_width" value="10*mm"/> + <constant name="Ecal_crystal_rotate_angle" value="30*degree"/> + <constant name="Ecal_crystal_envelope_length" value="266*mm"/> <!--not necessary, Ecal_barrel_inner_radius, Ecal_barrel_outer_radius and gap define the range of crystal--> + <constant name="numberZ" value="(int)Ecal_barrel_half_length*2/Ecal_crystal_y_width"/> + <constant name="Ecal_barrel_half_length_correct" value="numberZ*Ecal_crystal_y_width/2"/> <!--Must be n*Ecal_Crystal_y_width! --> + <constant name="Ecal_barrel_outer_radius_redef" value="sqrt(Ecal_barrel_inner_radius*Ecal_barrel_inner_radius+Ecal_crystal_envelope_length*Ecal_crystal_envelope_length + -2*Ecal_barrel_inner_radius*Ecal_crystal_envelope_length*cos(pi-Ecal_crystal_rotate_angle))"/> + </define> + + <regions> + <region name="EcalBarrelRegion"> + </region> + </regions> + + <detectors> + <detector id="DetID_ECAL" name="EcalBarrel" type="DD4hep_RotatedCrystalCalorimeter_v01" readout="EcalBarrelCollection" vis="Invisible" sensitive="true"> + <envelope> + <!--shape type="BooleanShape" operation="Union" material="Air"--> + <shape type="Tube" rmin="Ecal_barrel_inner_radius" rmax="Ecal_barrel_outer_radius_redef" dz="Ecal_barrel_half_length" material="Air"/> + <!--/shape--> + </envelope> + <dimensions rmin="Ecal_barrel_inner_radius" rmax="Ecal_barrel_outer_radius_redef" zhalf="Ecal_barrel_half_length_correct" alpha="Ecal_crystal_rotate_angle" nphi="1368" nz="numberZ" gap="0"/> + <crystal material="G4_BGO" vis="ECALVis"/> + </detector> + </detectors> + + <readouts> + <readout name="EcalBarrelCollection"> + <id>system:5,module:11,crystal:16</id> + </readout> + </readouts> + +</lccdd> diff --git a/Detector/DetCRD/compact/Standalone/Dimensions_v01_01.xml b/Detector/DetCRD/compact/Standalone/Dimensions_v01_01.xml new file mode 100644 index 0000000000000000000000000000000000000000..5b2a190c38b72997a83ad6377c1da2e017b07323 --- /dev/null +++ b/Detector/DetCRD/compact/Standalone/Dimensions_v01_01.xml @@ -0,0 +1,260 @@ +<?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="StandaloneDimensions" + title="master file with includes and world dimension" + author="C.D.Fu" + url="no" + status="development" + version="1.0"> + <comment> + undeterminded parameters + </comment> + </info> + + <define> + <constant name="CrossingAngle" value="0.033*rad"/> + + <constant name="Global_endcap_costheta" value="0.99"/> + + <constant name="GlobalTrackerReadoutID_DCH" type="string" value="system:8,chamber:1,layer:7,phi:16"/> + <constant name="GlobalTrackerReadoutID" type="string" value="system:5,side:-2,layer:9,module:8,sensor:8,barrelside:-2"/> + + <constant name="Field_nominal_value" value="3*tesla"/> + <constant name="Field_outer_nominal_value" value="-1.3*tesla"/> + + <constant name="env_safety" value="0.1*mm"/> + + <constant name="DetID_NOTUSED" value=" 0"/> + <constant name="DetID_VXD" value=" 1"/> + <constant name="DetID_SIT" value=" 2"/> + <constant name="DetID_FTD" value=" 3"/> + <constant name="DetID_TPC" value=" 4"/> + <constant name="DetID_SET" value=" 5"/> + <constant name="DetID_ETD" value=" 6"/> + <constant name="DetID_DC" value=" 4"/> <!--in order to cheat Clupatra, same as TPC--> + + <constant name="DetID_ECAL" value=" 20"/> + <constant name="DetID_ECAL_PLUG" value=" 21"/> + <constant name="DetID_HCAL" value=" 22"/> + <constant name="DetID_HCAL_RING" value=" 23"/> + <constant name="DetID_LCAL" value=" 24"/> + <constant name="DetID_BCAL" value=" 25"/> + <constant name="DetID_LHCAL" value=" 26"/> + <constant name="DetID_YOKE" value=" 27"/> + <constant name="DetID_COIL" value=" 28"/> + <constant name="DetID_ECAL_ENDCAP" value=" 29"/> + <constant name="DetID_HCAL_ENDCAP" value=" 30"/> + <constant name="DetID_YOKE_ENDCAP" value=" 31"/> + + <constant name="DetID_bwd" value="-1"/> + <constant name="DetID_barrel" value=" 0"/> + <constant name="DetID_fwd" value="+1"/> + + <constant name="BeamPipe_Be_inner_thickness" value="0.5*mm"/> + <constant name="BeamPipe_Cooling_thickness" value="0.5*mm"/> + <constant name="BeamPipe_Be_outer_thickness" value="0.3*mm"/> + <constant name="BeamPipe_Be_total_thickness" value="BeamPipe_Be_inner_thickness+BeamPipe_Cooling_thickness+BeamPipe_Be_outer_thickness"/> + <constant name="BeamPipe_Al_thickness" value="BeamPipe_Be_total_thickness"/> + <constant name="BeamPipe_Cu_thickness" value="2.0*mm"/> + + <constant name="BeamPipe_CentralBe_zmax" value="120*mm"/> + <constant name="BeamPipe_CentralAl_zmax" value="205*mm"/> + <constant name="BeamPipe_ConeAl_zmax" value="655*mm"/> + <constant name="BeamPipe_LinkerAl_zmax" value="700*mm"/> + <constant name="BeamPipe_LinkerCu_zmax" value="780*mm"/> + <constant name="BeamPipe_Waist_zmax" value="805*mm"/> + <constant name="BeamPipe_Crotch_zmax" value="855*mm"/> + <constant name="BeamPipe_FirstSeparated_zmax" value="1110*mm"/> + <constant name="BeamPipe_SecondSeparated_zmax" value="2200*mm"/> + <constant name="BeamPipe_end_z" value="12*m"/> + + <constant name="BeamPipe_Central_inner_radius" value="14*mm"/> + <constant name="BeamPipe_Expanded_inner_radius" value="20*mm"/> + <constant name="BeamPipe_Upstream_inner_radius" value="6*mm"/> + <constant name="BeamPipe_Dnstream_inner_radius" value="10*mm"/> + <constant name="BeamPipe_Crotch_hole_height" value="30.67*mm"/> + <constant name="BeamPipe_VertexRegion_rmax" value="BeamPipe_Central_inner_radius+BeamPipe_Al_thickness"/> + <constant name="BeamPipe_ForwardRegion_rmax" value="BeamPipe_Expanded_inner_radius+BeamPipe_Cu_thickness"/> + + <constant name="Vertex_inner_radius" value="BeamPipe_Central_inner_radius+BeamPipe_Be_total_thickness"/> + <constant name="Vertex_outer_radius" value="101*mm"/> + <constant name="Vertex_half_length" value="200*mm"/> + + <constant name="DC_Endcap_dz" value="0.1*mm"/> + <constant name="DC_half_length" value="2980*mm" /> + <constant name="DC_safe_distance" value="0.02*mm"/> + <constant name="SDT_inner_wall_thickness" value="0.2*mm"/> + <constant name="SDT_outer_wall_thickness" value="2.8*mm"/> + <constant name="MainTracker_half_length" value="DC_half_length+DC_Endcap_dz" /> + + <!--obselete for single drift chamber--> + <constant name="InnerTracker_half_length" value="DC_half_length" /> + <constant name="InnerTracker_inner_radius" value="234*mm"/> + <constant name="InnerTracker_outer_radius" value="909*mm"/> + <constant name="OuterTracker_half_length" value="DC_half_length"/> + <constant name="OuterTracker_inner_radius" value="1082.18*mm"/> + <constant name="OuterTracker_outer_radius" value="1723*mm"/> + + <!-- Parameters of single drift chamber --> + <constant name="DC_chamber_layer_rbegin" value="800*mm"/> + <constant name="DC_chamber_layer_rend" value="1800*mm"/> + + <constant name="DC_inner_radius" value="DC_chamber_layer_rbegin-SDT_inner_wall_thickness-DC_safe_distance"/> + <constant name="DC_outer_radius" value="DC_chamber_layer_rend+SDT_outer_wall_thickness+DC_safe_distance"/> + + <constant name="SIT1_inner_radius" value="230*mm"/> + <constant name="SIT2_inner_radius" value="410*mm"/> + <constant name="SIT3_inner_radius" value="590*mm"/> + <constant name="SIT4_inner_radius" value="770*mm"/> + <constant name="SIT1_half_length" value="461*mm"/> + <constant name="SIT2_half_length" value="691*mm"/> + <constant name="SIT3_half_length" value="1013*mm"/> + <constant name="SIT4_half_length" value="1335*mm"/> + + <constant name="SET_inner_radius" value="1815*mm"/> + + <constant name="SiTracker_barrel_endcap_gap" value="5*mm"/> + <constant name="SiTracker_DC_endcap_gap" value="10*mm"/> + <constant name="SiTracker_endcap_z1" value="SIT1_half_length+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_z2" value="SIT2_half_length+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_z3" value="SIT3_half_length+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_z4" value="SIT4_half_length+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_z5" value="MainTracker_half_length+SiTracker_DC_endcap_gap"/> + <constant name="SiTracker_endcap_outer_radius1" value="SIT1_inner_radius+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_outer_radius2" value="SIT2_inner_radius+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_outer_radius3" value="SIT3_inner_radius+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_outer_radius4" value="SIT4_inner_radius+SiTracker_barrel_endcap_gap"/> + <constant name="SiTracker_endcap_outer_radius5" value="SET_inner_radius+SiTracker_barrel_endcap_gap"/> + <!--obseleted --> + <constant name="FTD_BeamPipe_cable_clearance" value="10*mm"/> + <constant name="FTD_BeamPipe_gap" value="15*mm"/> + <constant name="FTD_InnerTracker_gap" value="5*mm"/> + + <!--obseleted constance, used by old construct, should be removed while creating new constrcut--> + <constant name="TPC_Ecal_Hcal_barrel_halfZ" value="MainTracker_half_length"/> + <constant name="TPC_inner_radius" value="InnerTracker_inner_radius"/> + <constant name="TPC_outer_radius" value="OuterTracker_outer_radius"/> + <constant name="SIT1_Radius" value="SIT1_inner_radius"/> + <constant name="SIT1_Half_Length_Z" value="SIT1_half_length"/> + <constant name="SIT2_Radius" value="InnerTracker_inner_radius"/> <!--fake, used by FTD_Simple_Staggered and FTD_cepc, now should be determined by inner tracker--> + <constant name="SIT2_Half_Length_Z" value="SIT2_half_length"/> + <constant name="TUBE_IPOuterTube_end_z" value="BeamPipe_CentralAl_zmax"/> + <constant name="TUBE_IPOuterTube_end_radius" value="BeamPipe_Central_inner_radius+BeamPipe_Al_thickness"/> + <constant name="TUBE_IPOuterBulge_end_z" value="BeamPipe_Crotch_zmax"/><!--"BeamPipe_ConeAl_zmax"/--> + <constant name="TUBE_IPOuterBulge_end_radius" value="BeamPipe_Crotch_zmax*tan(CrossingAngle/2)+BeamPipe_Dnstream_inner_radius+BeamPipe_Cu_thickness"/> + <!--"BeamPipe_Expanded_inner_radius+BeamPipe_Al_thickness+5*mm"/--> + + <constant name="Ecal_barrel_inner_radius" value="1900*mm"/> + <constant name="Ecal_barrel_thickness" value="280*mm"/> + <constant name="Ecal_barrel_outer_radius" value="(Ecal_barrel_inner_radius+Ecal_barrel_thickness)/cos(pi/8)"/> + <constant name="Ecal_barrel_half_length" value="3350*mm"/> + <constant name="Ecal_barrel_symmetry" value="8"/> + + <constant name="Ecal_endcap_inner_radius" value="350*mm"/> + <constant name="Ecal_endcap_outer_radius" value="Ecal_barrel_inner_radius+Ecal_barrel_thickness"/> + <constant name="Ecal_endcap_zmin" value="3050*mm"/> + <constant name="Ecal_endcap_zmax" value="3350*mm"/> + <constant name="Ecal_endcap_symmetry" value="8"/> + <!--obseleted constance, used by old construct, should be removed while creating new constrcut--> + <constant name="EcalEndcap_outer_radius" value="Ecal_barrel_outer_radius"/> + + <constant name="Solenoid_inner_radius" value="2330*mm"/> + <constant name="Solenoid_outer_radius" value="2480*mm"/> + <constant name="Solenoid_half_length" value="3830*mm"/> + <constant name="SolenoidCoil_half_length" value="3800*mm"/> + <constant name="SolenoidCoil_radius" value="2351*mm"/> + <constant name="SolenoidCoil_center_radius" value="(Solenoid_inner_radius+Solenoid_outer_radius)/2"/> + + <constant name="Hcal_barrel_inner_radius" value="2530*mm"/> + <constant name="Hcal_barrel_outer_radius" value="3610*mm"/> + <constant name="Hcal_barrel_half_length" value="4480*mm"/> + <constant name="Hcal_barrel_symmetry" value="12"/> + + <constant name="Hcal_endcap_inner_radius" value="400*mm"/> + <constant name="Hcal_endcap_outer_radius" value="Ecal_barrel_outer_radius"/> + <constant name="Hcal_endcap_zmin" value="3400*mm"/> + <constant name="Hcal_endcap_zmax" value="4480*mm"/> + <constant name="Hcal_endcap_symmetry" value="12"/> + <!--obseleted constance, used by old construct, should be removed while creating new constrcut--> + <constant name="HcalEndcap_max_z" value="Hcal_endcap_zmax"/> + <constant name="Hcal_endcap_outer_symmetry" value="Hcal_endcap_symmetry"/> + <constant name="Hcal_outer_radius" value="Hcal_endcap_outer_radius"/> + + <constant name="Yoke_barrel_inner_radius" value="3660*mm"/> + <constant name="Yoke_barrel_outer_radius" value="4260*mm"/> + <constant name="Yoke_barrel_half_length" value="Hcal_endcap_zmax"/> + <constant name="Yoke_barrel_symmetry" value="12"/> + + <constant name="Yoke_endcap_inner_radius" value="400*mm"/> + <constant name="Yoke_endcap_outer_radius" value="Yoke_barrel_outer_radius"/> + <constant name="Yoke_endcap_zmin" value="4660*mm"/> + <constant name="Yoke_endcap_zmax" value="5460*mm"/> + <constant name="Yoke_endcap_outer_symmetry" value="Yoke_barrel_symmetry"/> + <constant name="Yoke_endcap_inner_symmetry" value="0"/> + <!--obseleted constance, used by old construct, should be removed while creating new constrcut--> + <constant name="Yoke_Z_start_endcaps" value="Yoke_endcap_zmin"/> + + <constant name="tracker_region_zmax" value="OuterTracker_half_length"/> + <constant name="tracker_region_rmax" value="OuterTracker_outer_radius"/> + + </define> + + <limits> + <limitset name="cal_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + <limitset name="dc_limits"> + <limit name="step_length_max" particles="*" value="10.0" unit="mm" /> + </limitset> + <limitset name="tracker_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + </limits> + + <regions> + <region name="BeampipeRegion"/> + <region name="VertexRegion"/> + <region name="ForwardRegion"/> + </regions> + + <display> + <vis name="VXDVis" alpha="0.1" r="0.1" g=".5" b=".5" showDaughters="true" visible="true"/> + <vis name="VXDLayerVis" alpha="1.0" r="0.1" g=".5" b=".5" showDaughters="true" visible="true"/> + <vis name="VXDSupportVis" alpha="1.0" r="0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="FTDVis" alpha="1.0" r="0.5" g="0.87" b="0.11" showDaughters="true" visible="true"/> + <vis name="FTDSupportVis" alpha="1.0" r="0.3" g="0.3" b="1.0" showDaughters="true" visible="true"/> + <vis name="FTDSensitiveVis" alpha="1.0" r="0.3" g="0.5" b="1.0" showDaughters="true" visible="true"/> + <vis name="DCVis" alpha="1.0" r="0.96" g="0.64" b="0.90" showDaughters="true" visible="true"/> + <vis name="DCLayerVis" alpha="1.0" r="0.96" g="0.64" b="0.90" showDaughters="false" visible="true"/> + <vis name="SITVis" alpha="0.0" r="0.54" g="0.59" b="0.93" showDaughters="true" visible="false"/> + <vis name="SITSupportVis" alpha="1.0" r="0.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="SITSensitiveVis" alpha="1.0" r="0.67" g="0.99" b="0.78" showDaughters="false" visible="true"/> + <vis name="SETVis" alpha="0.0" r="0.8" g="0.8" b="0.4" showDaughters="true" visible="false"/> + <vis name="SETSupportVis" alpha="1.0" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="SETSensitiveVis" alpha="1.0" r="0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="ECALVis" alpha="1.0" r="0.2" g="0.6" b="0" showDaughters="true" visible="true"/> + <vis name="HCALVis" alpha="1.0" r="0.95" g="0.78" b="0.69" showDaughters="true" visible="true"/> + <vis name="SOLVis" alpha="1.0" r="0.4" g="0.4" b="0.4" showDaughters="true" visible="true"/> + <vis name="YOKEVis" alpha="1.0" r="0.64" g="0.75" b="0.99" showDaughters="false" visible="true"/> + <vis name="LCALVis" alpha="1.0" r="0.25" g="0.88" b="0.81" showDaughters="true" visible="true"/> + <vis name="SupportVis" alpha="1.0" r="0.2" g="0.2" b="0.2" showDaughters="true" visible="true"/> + <vis name="ShellVis" alpha="1.0" r="0.83" g="0.55" b="0.89" showDaughters="false" visible="true"/> + + <vis name="WhiteVis" alpha="0.0" r=".96" g=".96" b=".96" showDaughters="true" visible="true"/> + <vis name="LightGrayVis" alpha="0.0" r=".75" g=".75" b=".75" showDaughters="true" visible="true"/> + <vis name="Invisible" alpha="0.0" r="0.0" g="0.0" b="0.0" showDaughters="false" visible="false"/> + <vis name="SeeThrough" alpha="0.0" r="0.0" g="0.0" b="0.0" showDaughters="true" visible="false"/> + <vis name="RedVis" alpha="1.0" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="GreenVis" alpha="1.0" r="0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="BlueVis" alpha="1.0" r="0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="CyanVis" alpha="1.0" r="0.0" g="1.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="MagentaVis" alpha="1.0" r="1.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="YellowVis" alpha="1.0" r="1.0" g="1.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="BlackVis" alpha="1.0" r="0.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="GrayVis" alpha="1.0" r="0.5" g="0.5" b="0.5" showDaughters="true" visible="true"/> + </display> + +</lccdd> diff --git a/Detector/DetCRD/compact/Standalone/Standalone-EcalRotCrystal.xml b/Detector/DetCRD/compact/Standalone/Standalone-EcalRotCrystal.xml new file mode 100644 index 0000000000000000000000000000000000000000..771de36ea41e8279ed0f38d7a51818c1bded9f19 --- /dev/null +++ b/Detector/DetCRD/compact/Standalone/Standalone-EcalRotCrystal.xml @@ -0,0 +1,49 @@ +<?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="StandaloneEcalRotCrystal" + title="CepC standalone calorimeter with rotated crystal" + author="C.D.Fu" + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC detector simulation models used for detector study </comment> + </info> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="../CRD_common_v01/materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="6*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="./Dimensions_v01_01.xml"/> + + <include ref="../CRD_common_v01/Ecal_Rotated_Crystal_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/Standalone-Sim-RotCrystal.py b/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py new file mode 100644 index 0000000000000000000000000000000000000000..a424a1b3af88ed4b9929cf87ade65891277ce763 --- /dev/null +++ b/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py @@ -0,0 +1,105 @@ +#!/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 = "Standalone/Standalone-EcalRotCrystal.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 = ["e-"] +#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 = "RCEcalSim00.root" +out.outputCommands = ["keep *"] + +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [genalg, detsimalg, out], + EvtSel = 'NONE', + EvtMax = 10, + ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc], + OutputLevel=INFO +) diff --git a/Detector/DetCRD/src/Calorimeter/RotatedCrystalCalorimeter_v01_geo.cpp b/Detector/DetCRD/src/Calorimeter/RotatedCrystalCalorimeter_v01_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ec6e9fa5c229c348edccc36089d8cb59e24b7c1d --- /dev/null +++ b/Detector/DetCRD/src/Calorimeter/RotatedCrystalCalorimeter_v01_geo.cpp @@ -0,0 +1,128 @@ +//========================================================================== +// RotatedCrystalCalorimeter implementation +//-------------------------------------------------------------------------- +// rmin : inner radius (minimum of 8 vertex of crystal) +// rmax : outer radius (if gap=0, maximum of 8 vertex of crystal) +// alpha : rotate angle +// nphi : number of crystals rotated in phi direction +// nz : number of crystals in z direction +// gap : space in tail of crystal for SiPM, electronics, etc. +// Author: FU Chengdong, IHEP +//========================================================================== + +#include "DD4hep/DetFactoryHelper.h" +#include "XML/Utilities.h" +#include "XML/Layering.h" + +using namespace std; +using namespace dd4hep; +using namespace dd4hep::detail; + +static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) { + std::cout << "This is the RotatedCrystalCalorimeter_v01:" << std::endl; + + xml_det_t x_det = e; + xml_dim_t dim = x_det.dimensions(); + string det_name = x_det.nameStr(); + Material air = description.air(); + double zhalf = dim.zhalf(); + double rmin = dim.rmin(); + double rmax = dim.rmax(); + double alpha = dim.alpha(); + int nphi = dim.nphi(); + int nz = dim.nz(); + double gap = dim.gap(); + std::cout << "rmin = " << rmin/mm << std::endl + << "rmax = " << rmax/mm << std::endl + << "zhalf = " << zhalf/mm << std::endl + << "alpha = " << alpha/degree << std::endl + << "nphi = " << nphi << std::endl + << "nz = " << nz << std::endl; + + DetElement cal(det_name, x_det.id()); + if(alpha>90*dd4hep::degree){ + std::cout << "alpha>90 degree! return null cal!" << std::endl; + return cal; + } + + // --- create an envelope volume and position it into the world --------------------- + Volume envelope = dd4hep::xml::createPlacedEnvelope( description, e, cal ); + dd4hep::xml::setDetectorTypeFlag( e, cal ) ; + if( description.buildType() == BUILD_ENVELOPE ) return cal; + envelope.setVisAttributes(description, "SeeThrough"); + + DetElement module0(cal, "module1", x_det.id()); + + double yhalf = zhalf/nz; + Tube moduleTube(rmin, rmax, yhalf); + Volume moduleVol("module", moduleTube, air); + moduleVol.setVisAttributes(description, "SeeThrough"); + // + // | __C + // | __/ | + // | __/ | + // |__/ | + // B|___________|__E=C' of another crystal + // A| D + // | OA=rmin,OC=OE=rmax + // | /_BAD=alpha + // |O + double dphi = 2*M_PI/nphi; + double angleEAO = M_PI - alpha; + double angleAEO = std::asin(sin(angleEAO)/rmax*rmin); + double AE = rmax/sin(angleEAO)*sin(alpha-angleAEO); + double angleOpen = dphi; + double angleBottom = (M_PI-angleOpen)/2; + double disNeighbor = rmin*sin(dphi/2)*2; + //double x1 = disNeighbor/sin(angleBottom)*sin(angleEAO-angleBottom); + double x1 = disNeighbor/sin(angleBottom)*sin(angleOpen+angleBottom-alpha); + double edge = AE - disNeighbor/sin(angleBottom)*sin(alpha); + double length = edge*cos(angleOpen/2) - gap; + double x2 = x1 + length*tan(angleOpen/2)*2; + double center2A = sqrt(x1*x1/4+length*length/4); + double angleCenterAO = angleEAO+(M_PI-angleBottom-atan(length/x1)); + double center2O = sqrt(rmin*rmin+x1*x1/4+length*length/4-2*rmin*center2A*cos(angleCenterAO)); + double phi0Center = -asin(sin(angleCenterAO)/center2O*center2A); + std::cout << "Single crystal's open angle = " << angleOpen/degree << std::endl; + std::cout << "AE = " << AE/mm << ", /_AOE = " << (alpha-angleAEO)/degree << std::endl; + std::cout << "OA = " << rmin/mm << ", /_AEO = " << angleAEO/degree << std::endl; + std::cout << "OE = " << rmax/mm << ", /_EAO = " << angleEAO/degree << std::endl; + std::cout << "x1 = " << x1/mm << " x2 = " << x2/mm << " y = " << yhalf*2/mm << " l = " << length/mm << std::endl; + + xml_comp_t x_crystal = x_det.child( _Unicode(crystal) ); + Material crystalMat = description.material(x_crystal.materialStr()); + string crystalVis = x_crystal.visStr(); + + Trapezoid crystalTrap(x1/2, x2/2, yhalf, yhalf, length/2); + Volume crystalVol("crystal_whole", crystalTrap, crystalMat); + crystalVol.setVisAttributes(description, crystalVis); + + sens.setType("calorimeter"); + // TODO: crystal pack + crystalVol.setSensitiveDetector(sens); + for (xml_coll_t xp(x_crystal, _U(slice)); xp; ++xp){ + // + xml_comp_t x_slice = xp; + } + + for(int crystal_id=1; crystal_id<=nphi; crystal_id++){ + double angleRot = -alpha + dphi/2 + (crystal_id-1)*dphi; + double phiCenter = phi0Center + (crystal_id-1)*dphi; + Transform3D trafo(RotationZYX(0, angleRot+M_PI/2, M_PI/2), Translation3D(center2O*cos(phiCenter), center2O*sin(phiCenter), 0)); + PlacedVolume pv = moduleVol.placeVolume(crystalVol, trafo); + pv.addPhysVolID("crystal", crystal_id); + DetElement crystal(module0, _toString(crystal_id,"crystal%d"), x_det.id()); + crystal.setPlacement(pv); + } + + for (int module_id = 1; module_id <= nz; module_id++){ + DetElement module = module_id > 1 ? module0.clone(_toString(module_id,"module%d")) : module0; + PlacedVolume pv = envelope.placeVolume(moduleVol, Position(0,0,-zhalf+(2*module_id-1)*yhalf)); + pv.addPhysVolID("module", module_id); + module.setPlacement(pv); + } + + return cal; +} + +DECLARE_DETELEMENT(DD4hep_RotatedCrystalCalorimeter_v01, create_detector) diff --git a/Detector/DetCRD/src/Calorimeter/RotatedPolyhedraBarrelCalorimeter_v01_geo.cpp b/Detector/DetCRD/src/Calorimeter/RotatedPolyhedraBarrelCalorimeter_v01_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..723766ef891efba33ddd4bc357a81de5ea73760d --- /dev/null +++ b/Detector/DetCRD/src/Calorimeter/RotatedPolyhedraBarrelCalorimeter_v01_geo.cpp @@ -0,0 +1,169 @@ +//========================================================================== +// RotatedPolyhedraBarrelCalorimeter_v01 implementation +//-------------------------------------------------------------------------- +// Author: FU Chengdong, IHEP +//========================================================================== + +#include "DD4hep/DetFactoryHelper.h" +#include "XML/Utilities.h" +#include "XML/Layering.h" + +using namespace std; +using namespace dd4hep; +using namespace dd4hep::detail; + +static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) { + std::cout << "This is the RotatedPolyhedraBarrelCalorimeter_v01:" << std::endl; + + xml_det_t x_det = e; + Layering layering(x_det); + xml_comp_t staves = x_det.staves(); + xml_dim_t dim = x_det.dimensions(); + string det_name = x_det.nameStr(); + Material air = description.air(); + double totalThickness = layering.totalThickness(); + double gap = dim.gap();//xml_dim_t(x_det).gap(); + int numSides = dim.numsides(); + double zhalf = dim.zhalf(); + double rmin = dim.rmin(); + double phi0 = dim.phi0(); + double zpos = dim.zpos(); + std::cout << "rmin = " << rmin << std::endl + << "rmax = " << rmin+totalThickness << std::endl + << "nstave = " << numSides << std::endl + << "phi0 = " << phi0 << std::endl + << "zoffset = " << zpos << std::endl; + + DetElement cal(det_name, x_det.id()); + + // --- create an envelope volume and position it into the world --------------------- + Volume envelope = dd4hep::xml::createPlacedEnvelope( description, e, cal ) ; + dd4hep::xml::setDetectorTypeFlag( e, cal ) ; + if( description.buildType() == BUILD_ENVELOPE ) return cal ; + envelope.setVisAttributes(description, "SeeThrough"); + + // Add the subdetector envelope to the structure. + DetElement stave0(cal, "stave1", x_det.id()); + Material matStave = air; + std::string visStave = "SeeThrough"; + if ( staves ) { + visStave = staves.visStr(); + matStave = description.material(staves.materialStr()); + } + + double innerAngle = 2*M_PI/numSides; + double halfInnerAngle = innerAngle/2; + double tan_inner = std::tan(halfInnerAngle)*2; + double innerFaceLen = rmin*tan_inner + totalThickness/std::sin(innerAngle); + double outerFaceLen = (rmin+totalThickness)*tan_inner - totalThickness/std::sin(innerAngle); + double staveThickness = totalThickness; + + Trapezoid staveTrdOuter(innerFaceLen/2, outerFaceLen/2, zhalf, zhalf, staveThickness/2); + Volume staveOuterVol("stave_outer", staveTrdOuter, air); + staveOuterVol.setVisAttributes(description, "SeeThrough"); + + Trapezoid staveTrdInner(innerFaceLen/2 - gap/2, outerFaceLen/2 - gap/2, zhalf, zhalf, staveThickness/2); + Volume staveInnerVol("stave_inner", staveTrdInner, matStave); + staveInnerVol.setVisAttributes(description, visStave); + + double layerOuterAngle = (M_PI - innerAngle)/2; + double layerInnerAngle = (M_PI/2 - layerOuterAngle); + double layerSideAngle = M_PI - layerOuterAngle*2; + double layer_pos_z = -staveThickness/2; + double layer_dim_x = innerFaceLen/2 - gap/2; + int layer_num = 1; + + for (xml_coll_t xc(x_det, _U(layer)); xc; ++xc) { + xml_comp_t x_layer = xc; + int repeat = x_layer.repeat(); // Get number of times to repeat this layer. + const Layer* lay = layering.layer(layer_num - 1); // Get the layer from the layering engine. + // Loop over repeats for this layer. + for (int j = 0; j < repeat; j++) { + string layer_name = _toString(layer_num, "layer%d"); + double layer_thickness = lay->thickness(); + DetElement layer(stave0, layer_name, layer_num); + //### layeringExtension->setLayer(layer_num, layer, layerNormal); + + layer_dim_x -= layer_thickness/std::tan(layerSideAngle); + // Layer position in Z within the stave. + layer_pos_z += layer_thickness/2; + // Layer box & volume + Volume layer_vol(layer_name, Box(layer_dim_x, zhalf, layer_thickness/2), air); + + // Create the slices (sublayers) within the layer. + double slice_pos_z = -layer_thickness/2; + int slice_number = 1; + for (xml_coll_t xk(x_layer, _U(slice)); xk; ++xk) { + xml_comp_t x_slice = xk; + string slice_name = _toString(slice_number, "slice%d"); + double slice_thickness = x_slice.thickness(); + Material slice_material = description.material(x_slice.materialStr()); + DetElement slice(layer, slice_name, slice_number); + + slice_pos_z += slice_thickness/2; + // Slice volume & box + Volume slice_vol(slice_name, Box(layer_dim_x, zhalf, slice_thickness/2), slice_material); + + if (x_slice.isSensitive()) { + sens.setType("calorimeter"); + slice_vol.setSensitiveDetector(sens); + } + // Set region, limitset, and vis. + slice_vol.setAttributes(description, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr()); + // slice PlacedVolume + PlacedVolume slice_phv = layer_vol.placeVolume(slice_vol, Position(0, 0, slice_pos_z)); + slice_phv.addPhysVolID("slice", slice_number); + + slice.setPlacement(slice_phv); + // Increment Z position for next slice. + slice_pos_z += slice_thickness / 2; + // Increment slice number. + ++slice_number; + } + // Set region, limitset, and vis. + layer_vol.setAttributes(description, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr()); + + // Layer physical volume. + PlacedVolume layer_phv = staveInnerVol.placeVolume(layer_vol, Position(0, 0, layer_pos_z)); + layer_phv.addPhysVolID("layer", layer_num); + layer.setPlacement(layer_phv); + + // Increment the layer X dimension. + //layer_dim_x += layer_thickness * std::tan(layerInnerAngle); // * 2; + // Increment the layer Z position. + layer_pos_z += layer_thickness / 2; + // Increment the layer number. + ++layer_num; + } + } + + // Add stave inner physical volume to outer stave volume. + PlacedVolume pv = staveOuterVol.placeVolume(staveInnerVol, dd4hep::Position(gap/2,0,0)); + if ( staves ) { + stave0.setVisAttributes(description, staves.visStr(), staveInnerVol); + stave0.setVisAttributes(description, staves.visStr(), staveOuterVol); + } + // Place the staves. + double innerRotation = innerAngle; + double offsetRotation = phi0 - M_PI; + double sectCenterRadius = rmin + totalThickness/2; + double offset = totalThickness/std::sin(innerAngle)/2; + double rotX = M_PI / 2; + + for (int istave = 1; istave <= numSides; istave++){ + DetElement stave = istave > 1 ? stave0.clone(_toString(istave,"stave%d")) : stave0; + double rotY = offsetRotation - (istave-1)*innerRotation; + double posX = sectCenterRadius*std::sin(rotY) + offset*std::cos(rotY); + double posY = -sectCenterRadius*std::cos(rotY) + offset*std::sin(rotY); + Transform3D trafo(RotationZYX(0, rotY, rotX), Translation3D(posX, posY, zpos)); + PlacedVolume pv = envelope.placeVolume(staveOuterVol, trafo); + pv.addPhysVolID("stave", istave); + pv.addPhysVolID("system", cal.id()); + pv.addPhysVolID("barrel", 0); + stave.setPlacement(pv); + } + + return cal; +} + +DECLARE_DETELEMENT(DD4hep_RotatedPolyhedraBarrelCalorimeter_v01, create_detector)