Skip to content
Snippets Groups Projects
Commit 90d66de9 authored by FU Chengdong's avatar FU Chengdong
Browse files

add rotated polyhedra and crystal calorimeter

parent c4926c86
No related branches found
No related tags found
No related merge requests found
<?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>
<?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>
<?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="Standalone"
title="CepC reference detctor with coil inside Hcal, pixel SIT/SET"
author="C.D.Fu"
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"/>
</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>
#!/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
)
//==========================================================================
// 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;
//Translation3D tran(center2O*cos(phiCenter), center2O*sin(phiCenter), 0);
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)
DECLARE_DEPRECATED_DETELEMENT(RotatedCrystalCalorimeter_v01, create_detector)
//==========================================================================
// 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");
//Volume motherVol = description.pickMotherVolume(sdet);
//PolyhedraRegular polyhedra(numSides, rmin, rmin + totalThickness, zhalf*2);
//Volume wholeVol(det_name+"_Whole", polyhedra, air);
//wholeVol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
// 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;
//#### LayeringExtensionImpl* layeringExtension = new LayeringExtensionImpl();
//#### Position layerNormal(0,0,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;//-innerRotation/2;
double sectCenterRadius = rmin + totalThickness/2;
double offset = totalThickness/std::sin(innerAngle)/2;
double rotX = M_PI / 2;
//double rotY = -offsetRotation;
//double posX = -sectCenterRadius * std::sin(rotY);
//double posY = sectCenterRadius * std::cos(rotY);
for (int istave = 1; istave <= numSides; istave++){
//for (int istave = 1; istave <= 3; 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);
// Not a valid volID: pv.addPhysVolID("stave", 0);
pv.addPhysVolID("stave", istave);
stave.setPlacement(pv);
//cal.add(stave);
}
//placeStaves(cal, stave, rmin, numSides, totalThickness, envelopeVol, innerAngle, staveOuterVol);
// Set envelope volume attributes.
//envelope.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
//double z_offset = dim.hasAttr(_U(z_offset)) ? dim.z_offset() : 0.0;
//Transform3D transform(RotationZ(M_PI / numSides), Translation3D(0, 0, z_offset));
//PlacedVolume pv_whole = envelope.placeVolume(wholeVol, transform);
//pv_whole.addPhysVolID("system", cal.id());
//pv_whole.addPhysVolID("barrel", 0);
//cal.setPlacement(pv_whole);
//#### cal.addExtension<SubdetectorExtension>(new SubdetectorExtensionImpl(cal));
//#### cal.addExtension<LayeringExtension>(layeringExtension);
return cal;
}
DECLARE_DETELEMENT(DD4hep_RotatedPolyhedraBarrelCalorimeter_v01, create_detector)
DECLARE_DEPRECATED_DETELEMENT(RotatedPolyhedraBarrelCalorimeter_v01, create_detector)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment