diff --git a/Detector/DetCRD/compact/CRD_common_v01/OTKBarrel_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/OTKBarrel_v01_01.xml index bbef416d608bd750de1363d29afe95d389fe20d7..7a1840cb7f538c4538f99a4679d2f06e70d0dc39 100644 --- a/Detector/DetCRD/compact/CRD_common_v01/OTKBarrel_v01_01.xml +++ b/Detector/DetCRD/compact/CRD_common_v01/OTKBarrel_v01_01.xml @@ -10,7 +10,7 @@ </info> <define> <constant name="OTKBarrel_total_length" value="2*OTKBarrel_half_length" /> - <constant name="OTKBarrel_ladder_total_thickness" value="15*mm" /> + <constant name="OTKBarrel_ladder_total_thickness" value="10*mm" /> <constant name="OTKBarrel_ladder_total_width" value="160*mm" /> <constant name="OTKBarrel_ladder_total_length" value="OTKBarrel_total_length" /> <constant name="OTKBarrel_ladder_support_thickness" value="1*mm" /> @@ -28,18 +28,17 @@ <constant name="OTKBarrel_pcb_length" value="10*mm" /> <constant name="OTKBarrel_pcb_zgap" value="1*mm" /> <constant name="OTKBarrel_asic_thickness" value="500*um" /> - <constant name="OTKBarrel_asic_width" value="130*mm" /> - <constant name="OTKBarrel_asic_length" value="5*mm" /> + <constant name="OTKBarrel_asic_width" value="80*mm" /> + <constant name="OTKBarrel_asic_length" value="10*mm" /> <constant name="OTKBarrel_asic_zgap" value="1*mm" /> </define> <detectors> - <detector id="DetID_OTKBarrel" limits="otk_limits" name="OTKBarrel" type="SiTracker_otkbarrel_v01" vis="OTKBarrelVis" readout="OTKBarrelCollection" insideTrackingVolume="true"> + <!--detector id="DetID_OTKBarrel" limits="otk_limits" name="OTKBarrel" type="SiTracker_otkbarrel_v01" vis="OTKBarrelVis" readout="OTKBarrelCollection" insideTrackingVolume="true"--> + <detector id="DetID_OTKBarrel" name="OTKBarrel" type="SiTracker_otkbarrel_v01" vis="OTKBarrelVis" readout="OTKBarrelCollection" insideTrackingVolume="true"> <envelope> - <shape type="BooleanShape" operation="Union" material="Air" > - <shape type="Tube" rmin="OTKBarrel_inner_radius" rmax="OTKBarrel_outer_radius" dz="OTKBarrel_half_length" /> - </shape> + <shape type="Tube" rmin="OTKBarrel_inner_radius" rmax="OTKBarrel_outer_radius" dz="OTKBarrel_half_length" material="Air"/> </envelope> <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_STRIP "/> @@ -85,7 +84,7 @@ <readouts> <readout name="OTKBarrelCollection"> - <id>system:5,side:-2,layer:9,module:8,active:8,sensor:8</id> + <id>system:5,side:-2,layer:9,module:8,sensor:8,active:8</id> </readout> </readouts> </lccdd> diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml index 4adde286dc53d79d3bfedf007434869fc81f9a5b..25016db7d570eece30e62ae5b1876fcdf54fb7cc 100644 --- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_Dimensions_v01_01.xml @@ -34,8 +34,10 @@ <constant name="DetID_TPC" value=" 4"/> <constant name="DetID_SET" value=" 5"/> <constant name="DetID_ETD" value=" 6"/> - <constant name="DetID_OTKBarrel" value=" 8"/> - <constant name="DetID_OTKEndCap" value=" 9"/> + <constant name="DetID_ITKBarrel" value=" 2"/> + <constant name="DetID_ITKEndCap" value=" 3"/> + <constant name="DetID_OTKBarrel" value=" 5"/> + <constant name="DetID_OTKEndCap" value=" 6"/> <constant name="DetID_ECAL" value=" 20"/> <constant name="DetID_ECAL_PLUG" value=" 21"/> @@ -117,12 +119,12 @@ <!-- Parameters of time of flight tracker --> <constant name="OTKBarrel_inner_radius" value="1800*mm"/> - <constant name="OTKBarrel_outer_radius" value="1860*mm"/> + <constant name="OTKBarrel_outer_radius" value="1830*mm"/> <constant name="OTKBarrel_half_length" value="2940*mm"/> <constant name="OTKBarrelLayer1_half_length" value="OTKBarrel_half_length"/> <constant name="OTKBarrelLayer2_half_length" value="OTKBarrel_half_length"/> - <constant name="OTKBarrel1_inner_radius" value="1825*mm"/> - <constant name="OTKBarrel2_inner_radius" value="1845*mm"/> + <constant name="OTKBarrel1_inner_radius" value="1805*mm"/> + <constant name="OTKBarrel2_inner_radius" value="1820*mm"/> <constant name="SET_inner_radius" value="1800*mm"/> diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-oldVersion.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-oldVersion.xml new file mode 100644 index 0000000000000000000000000000000000000000..9578b78ec77543564994a7de57df2300d6c9adea --- /dev/null +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-oldVersion.xml @@ -0,0 +1,76 @@ +<?xml version="1.0" encoding="UTF-8"?> +<!-- option for old version of TDR_o1_v01 before frozen, to obsolete after testing new geometry--> +<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="TDR_o1_v01" + title="CepC reference detctor for TDR" + author="" + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC reference detector simulation models used for TDR </comment> + </info> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="../CRD_common_v02/materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="10*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="./TDR_Dimensions_v01_01.xml"/> + + <!--TODO: vertex cooling--> + <include ref="../CRD_common_v02/Beampipe_v01_03.xml"/> + <!--preliminary vertex and tracker, to update/--> + <include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/> + <include ref="../CRD_common_v02/FTD_SkewRing_v01_07.xml"/> + <include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/> + <!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/--> + <!--use 10 rows clustering version/--> + <include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/> + <include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/> + + <include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_02.xml"/> + <include ref="../CRD_common_v01/Ecal_Crystal_Endcap_v01_01.xml"/> + <include ref="../CRD_common_v01/SHcalGlass_Barrel_v05.xml"/> + <include ref="../CRD_common_v01/SHcalGlass_Endcaps_v01.xml"/> + + <!--Lumical to update--> + <include ref="../CRD_common_v01/Lumical_o1_v01.xml"/> + <!--preliminary Magnet, to update/--> + <include ref="../CRD_common_v02/Coil_Simple_v01_02.xml"/> + <!--preliminary Muon, obselete/--> + <!--include ref="../CRD_common_v02/Yoke_Polyhedra_Barrel_v01_01.xml"/> + <include ref="../CRD_common_v02/Yoke_Polyhedra_Endcaps_v01_01.xml"/--> + + <!--muon detector--> + <include ref="../CRD_common_v01/Muon_Barrel_v01_01.xml"/> + <include ref="../CRD_common_v01/Muon_Endcap_v01_02.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/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml index 83fba0bb1be521bca8794d53203a4b7e6a01535c..6dd00184a8e9c0a10e303c97def98944ad765849 100644 --- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml @@ -29,13 +29,16 @@ <!--TODO: vertex cooling--> <include ref="../CRD_common_v02/Beampipe_v01_03.xml"/> <!--preliminary vertex and tracker, to update/--> - <include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/> + <!--include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/--> + <include ref="../CRD_common_v02/VXD_Composite_v01_01.xml"/> <include ref="../CRD_common_v02/FTD_SkewRing_v01_07.xml"/> - <include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/> + <!--include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/--> + <include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/> <!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/--> <!-- use 10 rows clustering version--> <include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/> - <include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/> + <!--include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/--> + <include ref="../CRD_common_v01/OTKBarrel_v01_01.xml"/> <fields> <field name="InnerSolenoid" type="solenoid" diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-patch01_ITKE8.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-patch01_ITKE8.xml index 63f44cf77e7ba63b9ef98d3836de50e7350deb75..a5e04cec5bbd1d366ba2506c5dcc16e55e472421 100644 --- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-patch01_ITKE8.xml +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-patch01_ITKE8.xml @@ -30,13 +30,16 @@ <!--TODO: vertex cooling--> <include ref="../CRD_common_v02/Beampipe_v01_03.xml"/> <!--preliminary vertex and tracker, to update/--> - <include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/> + <!--include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/--> + <include ref="../CRD_common_v02/VXD_Composite_v01_01.xml"/> <include ref="../CRD_common_v02/FTD_SkewRing_v01_08.xml"/> - <include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/> + <!--include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/--> + <include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/> <!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/--> <!--use 10 rows clustering version/--> <include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/> - <include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/> + <!--include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/--> + <include ref="../CRD_common_v01/OTKBarrel_v01_01.xml"/> <include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_02.xml"/> <!--preliminary EcalEndcaps/--> diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml index 1cd781b9fd3344d91851b3ea24f63548f1b159fb..88dc0c6fe13bc9f83c050cdef76650fcdafd6e98 100644 --- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01.xml @@ -30,13 +30,16 @@ <!--TODO: vertex cooling--> <include ref="../CRD_common_v02/Beampipe_v01_03.xml"/> <!--preliminary vertex and tracker, to update/--> - <include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/> + <!--include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/--> + <include ref="../CRD_common_v02/VXD_Composite_v01_01.xml"/> <include ref="../CRD_common_v02/FTD_SkewRing_v01_07.xml"/> - <include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/> + <!--include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/--> + <include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/> <!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/--> <!--use 10 rows clustering version/--> <include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/> - <include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/> + <!--include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/--> + <include ref="../CRD_common_v01/OTKBarrel_v01_01.xml"/> <include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_02.xml"/> <include ref="../CRD_common_v01/Ecal_Crystal_Endcap_v01_01.xml"/> diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py index dbf60e86d54544ca9a44748c83582789032fc5c9..acf8c8a07cd76728a6bd1678b9630a3d0f30c97b 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py @@ -16,6 +16,7 @@ rndmgensvc = RndmGenSvc("RndmGenSvc") rndmgensvc.Engine = rndmengine.name() # option for standalone tracker study +#geometry_option = "TDR_o1_v01/TDR_o1_v01-oldVersion.xml" geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" if not os.getenv("DETCRDROOT"): diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking-oldVersion.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking-oldVersion.py new file mode 100644 index 0000000000000000000000000000000000000000..884cd89d82148db7bc81e867a7fc4b4ed2c16a2e --- /dev/null +++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking-oldVersion.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python +# option for old version of TDR_o1_v01 before frozen, to obsolete after testing new geometry +import os +from Gaudi.Configuration import * + +from Configurables import k4DataSvc +dsvc = k4DataSvc("EventDataSvc", input="sim00.root") + +from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ +seed = [12340] +# 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 = "TDR_o1_v01/TDR_o1_v01.xml-oldVersion" + +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 + +from Configurables import MarlinEvtSeeder +evtseeder = MarlinEvtSeeder("EventSeeder") + +from Configurables import GearSvc +gearsvc = GearSvc("GearSvc") + +from Configurables import TrackSystemSvc +tracksystemsvc = TrackSystemSvc("TrackSystemSvc") + +from Configurables import SimplePIDSvc +pidsvc = SimplePIDSvc("SimplePIDSvc") +cepcswdatatop = "/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest" +pidsvc.ParFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Service/SimplePIDSvc/data/dNdx_TPC.root") + +from Configurables import PodioInput +podioinput = PodioInput("PodioReader", collections=[ +# "EventHeader", + "MCParticle", + "VXDCollection", + "SITCollection", + "TPCCollection", + "SETCollection", + "FTDCollection" + ]) + +# digitization +vxdhitname = "VXDTrackerHits" +sithitname = "SITTrackerHits" +gashitname = "TPCTrackerHits" +sethitname = "SETTrackerHits" +setspname = "SETSpacePoints" +ftdhitname = "FTDTrackerHits" +ftdspname = "FTDSpacePoints" +from Configurables import PlanarDigiAlg +digiVXD = PlanarDigiAlg("VXDDigi") +digiVXD.SimTrackHitCollection = "VXDCollection" +digiVXD.TrackerHitCollection = vxdhitname +digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation" +digiVXD.ResolutionU = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004] +digiVXD.ResolutionV = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004] +digiVXD.UsePlanarTag = True +digiVXD.ParameterizeResolution = False +digiVXD.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +digiVXD.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +#digiVXD.OutputLevel = DEBUG + +digiSIT = PlanarDigiAlg("SITDigi") +digiSIT.IsStrip = False +digiSIT.SimTrackHitCollection = "SITCollection" +digiSIT.TrackerHitCollection = sithitname +digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation" +digiSIT.ResolutionU = [0.0098] +digiSIT.ResolutionV = [0.0433] +digiSIT.UsePlanarTag = True +digiSIT.ParameterizeResolution = False +digiSIT.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiSIT.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] +#digiSIT.OutputLevel = DEBUG + +digiSET = PlanarDigiAlg("SETDigi") +digiSET.IsStrip = False +digiSET.SimTrackHitCollection = "SETCollection" +digiSET.TrackerHitCollection = sethitname +digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation" +digiSET.ResolutionU = [0.005] +digiSET.ResolutionV = [0.021] +digiSET.UsePlanarTag = True +digiSET.ParameterizeResolution = False +digiSET.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiSET.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] +#digiSET.OutputLevel = DEBUG + +digiFTD = PlanarDigiAlg("FTDDigi") +digiFTD.IsStrip = False +digiFTD.SimTrackHitCollection = "FTDCollection" +digiFTD.TrackerHitCollection = ftdhitname +digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" +digiFTD.ResolutionU = [0.0072] +digiFTD.ResolutionV = [0.086] +digiFTD.UsePlanarTag = True +digiFTD.ParameterizeResolution = False +digiFTD.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiFTD.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] +#digiFTD.OutputLevel = DEBUG + +from Configurables import TPCDigiAlg +digiTPC = TPCDigiAlg("TPCDigi") +digiTPC.TPCCollection = "TPCCollection" +digiTPC.TPCLowPtCollection = "TPCLowPtCollection" +digiTPC.TPCTrackerHitsCol = gashitname +#default value, modify them according to future Garfield simulation results +#digiTPC.PixelClustering = True +#digiTPC.PointResolutionRPhi = 0.144 +#digiTPC.DiffusionCoeffRPhi = 0.0323 +#digiTPC.PointResolutionZ = 0.4 +#digiTPC.DiffusionCoeffZ = 0.23 +#digiTPC.N_eff = 30 +#digiTPC.OutputLevel = DEBUG + +# tracking +from Configurables import KalTestTool +# Close multiple scattering and smooth, used by clupatra +kt010 = KalTestTool("KalTest010") +kt010.MSOn = False +kt010.Smooth = False +#kt010.OutputLevel = DEBUG + +# Open multiple scattering, energy loss and smooth (default) +kt111 = KalTestTool("KalTest111") +#kt111.OutputLevel = DEBUG + +# Close smooth +kt110 = KalTestTool("KalTest110") +kt110.Smooth = False +#kt110.OutputLevel = DEBUG + +from Configurables import SiliconTrackingAlg +tracking = SiliconTrackingAlg("SiliconTracking") +tracking.LayerCombinationsFTD = [] +tracking.HeaderCol = "EventHeader" +tracking.VTXHitCollection = vxdhitname +tracking.SITHitCollection = sithitname +tracking.FTDPixelHitCollection = ftdhitname +tracking.FTDSpacePointCollection = ftdspname +tracking.SITRawHitCollection = "NotNeedForPixelSIT" +tracking.FTDRawHitCollection = ftdhitname +tracking.UseSIT = True +tracking.SmoothOn = False +tracking.NDivisionsInTheta = 10 +tracking.NDivisionsInPhi = 60 +tracking.NDivisionsInPhiFTD = 16 +tracking.MinDistCutAttach = 50 +# for p=1GeV, theta=10degree, Chi2FitCut = 1500, HelixMaxChi2 = 1000000, Chi2WZ = 0.02 +tracking.Chi2FitCut = 200 +tracking.MaxChi2PerHit = 200 +tracking.HelixMaxChi2 = 50000 +tracking.Chi2WZTriplet = 0.1 +tracking.Chi2WZQuartet = 0.1 +tracking.Chi2WZSeptet = 0.1 +#tracking.FitterTool = "KalTestTool/KalTest111" +#tracking.OutputLevel = DEBUG + +from Configurables import ForwardTrackingAlg +forward = ForwardTrackingAlg("ForwardTracking") +forward.FTDPixelHitCollection = ftdhitname +forward.FTDSpacePointCollection = ftdspname +forward.FTDRawHitCollection = ftdhitname +forward.Chi2ProbCut = 0.0 +forward.HitsPerTrackMin = 3 +forward.BestSubsetFinder = "SubsetSimple" +forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation", + "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"] +forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0] +forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05] +#forward.FitterTool = "KalTestTool/KalTest110" +#forward.OutputLevel = DEBUG + +from Configurables import TrackSubsetAlg +subset = TrackSubsetAlg("TrackSubset") +subset.TrackInputCollections = ["ForwardTracks", "SiTracks"] +subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, ftdspname] +subset.TrackSubsetCollection = "SubsetTracks" +#subset.FitterTool = "KalTestTool/KalTest111" +#subset.OutputLevel = DEBUG + +from Configurables import ClupatraAlg +clupatra = ClupatraAlg("Clupatra") +clupatra.TPCHitCollection = gashitname +#clupatra.OutputLevel = DEBUG + +from Configurables import FullLDCTrackingAlg +full = FullLDCTrackingAlg("FullTracking") +full.VTXTrackerHits = vxdhitname +full.SITTrackerHits = sithitname +full.TPCTrackerHits = gashitname +full.SETTrackerHits = sethitname +full.FTDPixelTrackerHits = ftdhitname +full.FTDSpacePoints = ftdspname +full.SITRawHits = "NotNeedForPixelSIT" +full.SETRawHits = "NotNeedForPixelSET" +full.FTDRawHits = ftdhitname +full.TPCTracks = "ClupatraTracks" # add standalone TPC track +full.SiTracks = "SubsetTracks" +full.OutputTracks = "CompleteTracks" # default name +full.VTXHitToTrackDistance = 5. +full.FTDHitToTrackDistance = 5. +full.SITHitToTrackDistance = 3. +full.SETHitToTrackDistance = 5. +full.MinChi2ProbForSiliconTracks = 0 +full.MaxChi2PerHit = 500 +#full.OutputLevel = DEBUG + +from Configurables import TPCDndxAlg +tpc_dndx = TPCDndxAlg("TPCDndxAlg") +tpc_dndx.Method = "Simple" + +from Configurables import TrackParticleRelationAlg +tpr = TrackParticleRelationAlg("Track2Particle") +tpr.MCParticleCollection = "MCParticle" +tpr.TrackList = ["CompleteTracks", "ClupatraTracks"] +tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] +#tpr.OutputLevel = DEBUG + +from Configurables import TrueMuonTagAlg +tmt = TrueMuonTagAlg("TrueMuonTag") +tmt.MCParticleCollection = "MCParticle" +tmt.TrackList = ["CompleteTracks"] +tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] +tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%) +tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta) +#tmt.OutputLevel = DEBUG + +# output +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "rec00.root" +out.outputCommands = ["keep *"] + +# ApplicationMgr +from Configurables import ApplicationMgr +mgr = ApplicationMgr( + TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tmt, out], + EvtSel = 'NONE', + EvtMax = 5, + ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc], + HistogramPersistency = 'ROOT', + OutputLevel = ERROR +) diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py index b621a182c161db2b19bc63dd6de62bbe2617f533..5c9b07cd8bce345471a6fa4f92c422410523b344 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py @@ -7,8 +7,8 @@ dsvc = k4DataSvc("EventDataSvc", input="sim00.root") from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ seed = [12340] -# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi -rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4 +# 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 @@ -51,7 +51,8 @@ podioinput = PodioInput("PodioReader", collections=[ "VXDCollection", "SITCollection", "TPCCollection", - "SETCollection", +# "SETCollection", + "OTKBarrelCollection", "FTDCollection" ]) @@ -59,23 +60,29 @@ podioinput = PodioInput("PodioReader", collections=[ vxdhitname = "VXDTrackerHits" sithitname = "SITTrackerHits" gashitname = "TPCTrackerHits" -sethitname = "SETTrackerHits" -setspname = "SETSpacePoints" +sethitname = "OTKBarrelTrackerHits" +setspname = "OTKBarrelSpacePoints" ftdhitname = "FTDTrackerHits" ftdspname = "FTDSpacePoints" -from Configurables import PlanarDigiAlg -digiVXD = PlanarDigiAlg("VXDDigi") +from Configurables import SmearDigiTool +vxdtool = SmearDigiTool("VXD") +vxdtool.ResolutionU = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004] +vxdtool.ResolutionV = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004] +vxdtool.UsePlanarTag = True +vxdtool.ParameterizeResolution = False +vxdtool.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +vxdtool.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +#vxdtool.OutputLevel = DEBUG + +from Configurables import SiTrackerDigiAlg +digiVXD = SiTrackerDigiAlg("VXDDigi") digiVXD.SimTrackHitCollection = "VXDCollection" digiVXD.TrackerHitCollection = vxdhitname digiVXD.TrackerHitAssociationCollection = "VXDTrackerHitAssociation" -digiVXD.ResolutionU = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004] -digiVXD.ResolutionV = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004] -digiVXD.UsePlanarTag = True -digiVXD.ParameterizeResolution = False -digiVXD.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] -digiVXD.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +digiVXD.DigiTool = "SmearDigiTool/VXD" #digiVXD.OutputLevel = DEBUG +from Configurables import PlanarDigiAlg digiSIT = PlanarDigiAlg("SITDigi") digiSIT.IsStrip = False digiSIT.SimTrackHitCollection = "SITCollection" @@ -91,9 +98,9 @@ digiSIT.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.296 digiSET = PlanarDigiAlg("SETDigi") digiSET.IsStrip = False -digiSET.SimTrackHitCollection = "SETCollection" +digiSET.SimTrackHitCollection = "OTKBarrelCollection" digiSET.TrackerHitCollection = sethitname -digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation" +digiSET.TrackerHitAssociationCollection = "OTKBarrelTrackerHitAssociation" digiSET.ResolutionU = [0.005] digiSET.ResolutionV = [0.021] digiSET.UsePlanarTag = True @@ -230,14 +237,14 @@ from Configurables import TrackParticleRelationAlg tpr = TrackParticleRelationAlg("Track2Particle") tpr.MCParticleCollection = "MCParticle" tpr.TrackList = ["CompleteTracks", "ClupatraTracks"] -tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] +tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] #tpr.OutputLevel = DEBUG from Configurables import TrueMuonTagAlg tmt = TrueMuonTagAlg("TrueMuonTag") tmt.MCParticleCollection = "MCParticle" tmt.TrackList = ["CompleteTracks"] -tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] +tmt.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "OTKBarrelTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"] tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%) tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta) #tmt.OutputLevel = DEBUG diff --git a/Detector/DetCRD/src/Tracker/SiTracker_itkbarrel_v02_geo.cpp b/Detector/DetCRD/src/Tracker/SiTracker_itkbarrel_v02_geo.cpp index 6706f477aa3a9668159d03decb3d10b4c77d7746..51bd66e56b97ee19ef0034ba79574fea5a5a2713 100644 --- a/Detector/DetCRD/src/Tracker/SiTracker_itkbarrel_v02_geo.cpp +++ b/Detector/DetCRD/src/Tracker/SiTracker_itkbarrel_v02_geo.cpp @@ -670,6 +670,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h pv = StaveLogical.placeVolume(TitubeLogical, Position(-stave_thickness/2.0 + tube_outer_radius, -support_width / 4.0, 0.)); pv = StaveLogical.placeVolume(TitubeLogical, Position(-stave_thickness/2.0 + tube_outer_radius, support_width / 4.0, 0.)); + double stave_radius = stave_sens_radius - stave_thickness/2 + (max_connector_thickness + flex_thickness + sensor_thickness/2); for(int i = 0; i < n_staves; i++){ std::stringstream stave_enum; stave_enum << "sit_stave_" << layer_id << "_" << i; @@ -706,7 +707,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h } } - double stave_radius = stave_sens_radius - stave_thickness/2 + (max_connector_thickness + flex_thickness + sensor_thickness/2); + //double stave_radius = stave_sens_radius - stave_thickness/2 + (max_connector_thickness + flex_thickness + sensor_thickness/2); Transform3D tr (RotationZYX(stave_dphi*i,0.,0.),Position(stave_radius*cos(stave_phi0+stave_dphi*i), stave_radius*sin(stave_phi0+stave_dphi*i), 0.)); //std::cout << "1st Test!!!" << endl; pv = layer_assembly.placeVolume(StaveLogical,tr); @@ -725,17 +726,19 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h Layer.ladderNumber = n_staves; Layer.phi0 = 0.; //Layer.modulesPerStave = n_modules_per_stave; - Layer.sensorsPerLadder = n_modules_per_stave*n_sensors_per_module; + Layer.sensorsPerLadder = n_modules_per_stave; + //Layer.sensorsPerLadder = n_modules_per_stave*n_sensors_per_module; //Layer.lengthModule = module_active_length; - Layer.lengthSensor = sensor_active_length; - Layer.distanceSupport = sensitive_radius; + Layer.lengthSensor = module_length;//sensor_active_length; + Layer.distanceSupport = stave_radius*cos(stave_phi0) + StaveSupportenv_start_height; //sensitive_radius; Layer.thicknessSupport = support_thickness / 2.0; - Layer.offsetSupport = -stave_offset; + Layer.offsetSupport = stave_radius*sin(stave_phi0);//-stave_offset; Layer.widthSupport = support_width; Layer.zHalfSupport = support_half_length; - Layer.distanceSensitive = sensitive_radius + tube_outer_radius*2. + support_thickness; - Layer.thicknessSensitive = module_thickness; - Layer.offsetSensitive = -stave_offset/2.0; + Layer.distanceSensitive = stave_radius*cos(stave_phi0) + StaveSupportenv_start_height + support_thickness; + //sensitive_radius + tube_outer_radius*2. + support_thickness; + Layer.thicknessSensitive = sensor_thickness;//module_thickness; + Layer.offsetSensitive = stave_radius*sin(stave_phi0);//-stave_offset;//stave_offset/2.0; Layer.widthSensitive = module_width; Layer.zHalfSensitive = support_half_length; diff --git a/Detector/DetCRD/src/Tracker/SiTracker_otkbarrel_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTracker_otkbarrel_v01_geo.cpp index b1f5ed127561eb6304c627a13ac5c58e58340433..45538fe223b2cafc9285cab8cc2b4c8f0cab3acc 100644 --- a/Detector/DetCRD/src/Tracker/SiTracker_otkbarrel_v01_geo.cpp +++ b/Detector/DetCRD/src/Tracker/SiTracker_otkbarrel_v01_geo.cpp @@ -71,7 +71,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h dd4hep::PlacedVolume pv; int layer_id = x_layer.attr<int>(_Unicode(layer_id)); - std::cout << "layer_id: " << layer_id << endl; + std::cout << "layer_id: " << layer_id << " ......." << endl; double sensitive_radius = x_layer.attr<double>(_Unicode(ladder_radius)); int n_ladders = x_layer.attr<int>(_Unicode(n_ladders)) ; @@ -233,7 +233,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h } //place the flex envelope inside the ladder envelope - pv = LadderLogical.placeVolume(FlexEnvelopeLogical, Position((support_height+flex_thickness)/2.0+sensor_thickness, (-support_width+flex_width)/2.0, 0.)); + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, Position((support_height+flex_thickness)/2.0+sensor_thickness+pcb_thickness+asic_thickness, (-support_width+flex_width)/2.0, 0.)); //create sensor envelope logical volume Box SensorEnvelopeSolid((sensor_thickness+pcb_thickness+asic_thickness) / 2.0, sensor_total_width / 2.0, support_length / 2.0); @@ -313,8 +313,8 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h pv = LadderSupportEnvelopeLogical.placeVolume(LadderSupportLogical); pv = LadderLogical.placeVolume(LadderSupportEnvelopeLogical); + float rot = layer_id*0.5; for(int i = 0; i < n_ladders; i++){ - float rot = layer_id*0.5; std::stringstream ladder_enum; ladder_enum << "otkbarrel_ladder_" << layer_id << "_" << i; DetElement ladderDE(layerDE, ladder_enum.str(), x_det.id()); @@ -347,7 +347,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h pv = layer_assembly.placeVolume(LadderLogical,tr); pv.addPhysVolID("module", i ) ; ladderDE.setPlacement(pv); - std::cout << ladder_enum.str() << " done." << endl; + //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; } @@ -355,19 +355,21 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h dd4hep::rec::ZPlanarData::LayerLayout otkbarrelLayer; otkbarrelLayer.ladderNumber = n_ladders; - otkbarrelLayer.phi0 = ladder_phi0; + otkbarrelLayer.phi0 = ladder_phi0 + ladder_dphi*rot; otkbarrelLayer.sensorsPerLadder = n_sensors_per_side; - otkbarrelLayer.lengthSensor = sensor_length; - otkbarrelLayer.distanceSupport = sensitive_radius; + otkbarrelLayer.lengthSensor = sensor_length + dead_gap; // dead region also has silicon material + otkbarrelLayer.distanceSupport = ladder_radius*cos(ladder_phi0) - support_thickness/2.0; //sensitive_radius; otkbarrelLayer.thicknessSupport = support_thickness / 2.0; otkbarrelLayer.offsetSupport = -ladder_offset; otkbarrelLayer.widthSupport = support_width; otkbarrelLayer.zHalfSupport = support_length / 2.0; - otkbarrelLayer.distanceSensitive = sensitive_radius + support_height / 2.0 + flex_thickness; + otkbarrelLayer.distanceSensitive = ladder_radius*cos(ladder_phi0) - support_thickness/2.0 + - sensor_thickness; //sensitive_radius + support_height / 2.0 + flex_thickness; otkbarrelLayer.thicknessSensitive = sensor_thickness; otkbarrelLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); otkbarrelLayer.widthSensitive = sensor_active_width; - otkbarrelLayer.zHalfSensitive = (n_sensors_per_side*(sensor_length + dead_gap) - dead_gap) / 2.0; + //otkbarrelLayer.zHalfSensitive = (n_sensors_per_side*(sensor_length + dead_gap) - dead_gap) / 2.0; + otkbarrelLayer.zHalfSensitive = n_sensors_per_side*(sensor_length + dead_gap) / 2.0; // add dead_gap to same sensor, little effect? zPlanarData->layers.push_back(otkbarrelLayer); } diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp index 5c1a9d9cd2bc924a1f15eb664d01efe5a75fc298..7b7ffd872a39fc7b890ac9085ea45a8b09563594 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp @@ -117,7 +117,7 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc declareProperty("VTXTrackerHits", _VTXTrackerHitColHdl, "VTX Hit Collection"); declareProperty("SITTrackerHits", _SITTrackerHitColHdl, "SIT Hit Collection"); declareProperty("SETTrackerHits", _SETTrackerHitColHdl, "SET Hit Collection"); - //declareProperty("ETDTrackerHits", _ETDTrackerHitColHdl, "ETD Hit Collection"); + declareProperty("ETDTrackerHits", _ETDTrackerHitColHdl, "ETD Hit Collection"); declareProperty("TPCTrackerHits", _TPCTrackerHitColHdl, "TPC Hit Collection"); declareProperty("SITRawHits", _inSITRawColHdl, "SIT Raw Hit Collection of SpacePoints"); declareProperty("SETRawHits", _inSETRawColHdl, "SET Raw Hit Collection of SpacePoints"); @@ -449,12 +449,16 @@ fitstart: trkHits.push_back(it->second); } - //for(int iHit=0;iHit<trkHits.size();iHit++){ - //const edm4hep::Vector3d& pos = trkHits[iHit].getPosition(); - //debug() << "Hit " << iHit << ": r= " << sqrt(pos.x*pos.x+pos.y*pos.y) << " id = " << trkHits[iHit].id() << endmsg; - //} - - debug() << "trkHits ready " << trkHits.size() << endmsg; + if (msgLevel(MSG::DEBUG)) { + debug() << "trkHits ready " << trkHits.size() << endmsg; + for (unsigned iHit=0; iHit<trkHits.size(); iHit++) { + const edm4hep::TrackerHit hit = trkHits[iHit]; + unsigned long cellID = hit.getCellID(); + const edm4hep::Vector3d& pos = hit.getPosition(); + double r = sqrt(pos.x*pos.x+pos.y*pos.y); + debug() << "Hit " << iHit << ": r= " << r << " id = " << hit.id() << " cellID = " << cellID << " det = " << (cellID&0x1F) << endmsg; + } + } int error_code = 0; edm4hep::MutableTrack track; @@ -1170,7 +1174,6 @@ void FullLDCTrackingAlg::prepareVectors() { _allSETHits.push_back( hitExt ); mapTrackerHits[trkhit] = hitExt; - double pos[3]; for (int i=0; i<3; ++i) { @@ -1180,7 +1183,139 @@ void FullLDCTrackingAlg::prepareVectors() { debug() << " SET Hit " << trkhit.id() << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " layer = " << layer << endmsg; } } - + + const edm4hep::TrackerHitCollection* hitETDCol = nullptr; + try { + hitETDCol = _ETDTrackerHitColHdl.get(); + } + catch ( GaudiException &e ) { + debug() << "Collection " << _ETDTrackerHitColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg; + } + + const edm4hep::TrackerHitCollection* rawETDCol = nullptr; + if(hitETDCol){ + try{ + rawETDCol = _inETDRawColHdl.get(); + } + catch ( GaudiException &e ) { + warning() << "Collection " << _inETDRawColHdl.fullKey() << " is unavailable in event " << _nEvt << ", if SIT is Space Point, it needed " << endmsg; + } + } + + if(hitETDCol){ + if(rawETDCol) Navigation::Instance()->AddTrackerHitCollection(rawETDCol); + + int nelem = hitETDCol->size(); + + debug() << "Number of ETD hits = " << nelem << endmsg; + + double drphi(NAN); + double dz(NAN); + + for(edm4hep::TrackerHit trkhit : *hitETDCol){ + // hit could be of the following type + // 1) TrackerHit, either ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT or just standard TrackerHit + // 2) TrackerHitPlane, either 1D or 2D + // 3) TrackerHitZCylinder, if coming from a simple cylinder design as in the LOI + + // Establish which of these it is in the following order of likelyhood + // i) ILDTrkHitTypeBit::ONE_DIMENSIONAL (TrackerHitPlane) Should Never Happen: SpacePoints Must be Used Instead + // ii) ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT (TrackerHit) + // iii) TrackerHitPlane (Two dimentional) + // iv) TrackerHitZCylinder + // v) Must be standard TrackerHit + int layer = getLayerID(trkhit); + + if (layer < 0 || (unsigned)layer >= _nLayersETD) { + fatal() << "FullLDCTrackingAlg => fatal error in ETD : layer is outside allowed range : " << layer << endmsg; + exit(1); + } + + // first check that we have not been given 1D hits by mistake, as they won't work here + if ( UTIL::BitSet32( trkhit.getType() )[ UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ] ) { + fatal() << "SiliconTrackingAlg => fatal error in SIT : layer is outside allowed range : " << layer << endmsg; + exit(1); + } + // most likely case: COMPOSITE_SPACEPOINT hits formed from stereo strip hits + else if ( UTIL::BitSet32( trkhit.getType() )[ UTIL::ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT ] ) { + // SJA:FIXME: fudge for now by a factor of two and ignore covariance + drphi = 2 * sqrt(trkhit.getCovMatrix()[0] + trkhit.getCovMatrix()[2]); + dz = sqrt(trkhit.getCovMatrix()[5]); + } + // or a PIXEL based SET, using 2D TrackerHitPlane like the VXD above + else if ( UTIL::BitSet32( trkhit.getType() )[3] ) { + // FIXME should consider this carefully + + // first we need to check if the measurement vectors are aligned with the global coordinates + // gear::Vector3D U(1.0,trkhit_P->getU()[1],trkhit_P->getU()[0],gear::Vector3D::spherical); + // gear::Vector3D V(1.0,trkhit_P->getV()[1],trkhit_P->getV()[0],gear::Vector3D::spherical); + // FIXME Should calculate it correctly + gear::Vector3D U(1.0,trkhit.getCovMatrix()[1],trkhit.getCovMatrix()[0],gear::Vector3D::spherical); + gear::Vector3D V(1.0,trkhit.getCovMatrix()[4],trkhit.getCovMatrix()[3],gear::Vector3D::spherical); + gear::Vector3D Z(0.0,0.0,1.0); + + const float eps = 1.0e-07; + // V must be the global z axis + if( fabs(1.0 - V.dot(Z)) > eps ) { + fatal() << "PIXEL ETD Hit measurment vectors V is not equal to the global Z axis. \n exit(1) called from file " << __FILE__ << " : " << __LINE__ << endmsg; + exit(1); + } + + // U must be normal to the global z axis + if( fabs(U.dot(Z)) > eps ) { + fatal() << "PIXEL ETD Hit measurment vectors U is not in the global X-Y plane. \n exit(1) called from file " << __FILE__ << " : " << __LINE__ << endmsg; + exit(1); + } + + // FIXME should use the correct + // drphi = trkhit_P->getdU(); + // dz = trkhit_P->getdV(); + drphi = trkhit.getCovMatrix()[2]; + dz = trkhit.getCovMatrix()[5]; + } + // or a simple cylindrical design, as used in the LOI + /*FIXME, fucd + else if ( true ) { + trkhit_C = hitCollection->at( ielem ); + // FIXME + // drphi = trkhit_C->getdRPhi(); + // dz = trkhit_C->getdZ(); + drphi = 1.0; + dz = 1.0; + } + */ + // this would be very unlikely, but who knows ... just an ordinary TrackerHit, which is not a COMPOSITE_SPACEPOINT + else { + // SJA:FIXME: fudge for now by a factor of two and ignore covariance + drphi = 2 * sqrt(trkhit.getCovMatrix()[0] + trkhit.getCovMatrix()[2]); + dz = sqrt(trkhit.getCovMatrix()[5]); + } + + // now that the hit type has been established carry on and create a + + TrackerHitExtended * hitExt = new TrackerHitExtended( trkhit ); + + // SJA:FIXME: just use planar res for now + hitExt->setResolutionRPhi(drphi); + hitExt->setResolutionZ(dz); + + // set type is now only used in one place where it is set to 0 to reject hits from a fit, set to INT_MAX to try and catch any missuse + hitExt->setType(int(INT_MAX)); + // det is no longer used set to INT_MAX to try and catch any missuse + hitExt->setDet(int(INT_MAX)); + + _allETDHits.push_back( hitExt ); + mapTrackerHits[trkhit] = hitExt; + + double pos[3]; + for (int i=0; i<3; ++i) { + pos[i] = trkhit.getPosition()[i]; + } + + debug() << " ETD Hit " << trkhit.id() << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " layer = " << layer << endmsg; + } + } + // Reading VTX Hits const edm4hep::TrackerHitCollection* hitVTXCol = nullptr; try { @@ -3348,94 +3483,56 @@ void FullLDCTrackingAlg::AddNotAssignedHits() { // currently any non-combined Silicon or TPC tracks are added to the list of tracks meaning their hits will not be available to be picked up. // it might be preferable to drop these tracks, at least for track in Silicon with very bad Chi2 probabilities, and allow their hits to be re-used here ... + if (_assignSETHits>0 || _assignETDHits>0) { + // Creating helix extrapolations + CreateExtrapolations(); + + if (_assignETDHits>0) { // Assignment of ETD Hits + debug() << "Assign ETD/OTKEndcap hits *********************************" << endmsg; + + int nETDHits = _allETDHits.size(); + std::vector<TrackerHitExtendedVec> ETDHits; + ETDHits.resize(_nLayersETD); + + for (int iETD=0;iETD<nETDHits;++iETD) { + TrackerHitExtended * trkHitExt = _allETDHits[iETD]; + edm4hep::TrackerHit trkHit = trkHitExt->getTrackerHit(); + int layer = getLayerID(trkHit); + if (layer >= 0 && layer < _nLayersETD) + ETDHits[layer].push_back(trkHitExt); + } + for (int iL=0; iL<_nLayersETD; ++iL) { // loop over ETD layers + TrackerHitExtendedVec hitVec = ETDHits[iL]; + int refit = 0; + if (hitVec.empty() == false) AssignOuterHitsToTracks(hitVec, _distCutForETDHits, refit); + } + } - - // Creating helix extrapolations - // if (_assignSETHits>0||_assignETDHits>0) - // CreateExtrapolations(); - - // if (_assignSETHits>0) { // Assignment of SET Hits - // - // const gear::GearParameters& pSETDet = Global::GEAR->getGearParameters("SET"); - // int nLayersSET = int(pSETDet.getDoubleVals("SETLayerRadius").size()); - // - // int nSETHits = _allSETHits.size(); - // std::vector<TrackerHitExtendedVec> SETHits; - // SETHits.resize(nLayersSET); - // - // for (int iSET=0;iSET<nSETHits;++iSET) { - // TrackerHitExtended * trkHit = _allSETHits[iSET]; - // edm4hep::TrackerHit hit = trkHit->getTrackerHit(); - // int layer = getLayerID(trkHit); - // if (layer>=0&&layer<nLayersSET) - // SETHits[layer].push_back(trkHit); - // } - // for (int iL=0; iL<nLayersSET; ++iL) { // loop over SET layers - // TrackerHitExtendedVec hitVec = SETHits[iL]; - // int refit = 1; - // AssignOuterHitsToTracks(hitVec,_distCutForSETHits,refit); - // } - // } - - // if (_assignETDHits>0) { // Assignment of ETD Hits - // - // const gear::GearParameters& pETDDet = Global::GEAR->getGearParameters("ETD"); - // int nLayersETD = int(pETDDet.getDoubleVals("ETDLayerZ").size()); - // - // int nETDHits = _allETDHits.size(); - // std::vector<TrackerHitExtendedVec> ETDHits; - // ETDHits.resize(nLayersETD); - // - // for (int iETD=0;iETD<nETDHits;++iETD) { - // TrackerHitExtended * trkHit = _allETDHits[iETD]; - // edm4hep::TrackerHit hit = trkHit->getTrackerHit(); - // int layer = getLayerID(trkHit); - // if (layer>=0 && layer < nLayersETD) - // ETDHits[layer].push_back(trkHit); - // } - // for (int iL=0; iL<nLayersETD; ++iL) { // loop over ETD layers - // TrackerHitExtendedVec hitVec = ETDHits[iL]; - // int refit = 0; - // AssignOuterHitsToTracks( hitVec, _distCutForETDHits, refit ); - // } - // - // } - - // // Cleaning up helix extrapolations - // if (_assignSETHits>0||_assignETDHits>0) - // CleanUpExtrapolations(); - - - - - - if (_assignSETHits>0) { // Assignment of SET Hits - debug() << "Assign SET hits *********************************" << endmsg; - - // Creating helix extrapolations - CreateExtrapolations(); + if (_assignSETHits>0) { // Assignment of SET Hits + debug() << "Assign SET hits *********************************" << endmsg; - int nSETHits = _allSETHits.size(); - std::vector<TrackerHitExtendedVec> SETHits; + int nSETHits = _allSETHits.size(); + std::vector<TrackerHitExtendedVec> SETHits; - SETHits.resize(_nLayersSET); + SETHits.resize(_nLayersSET); - for (int iSET=0;iSET<nSETHits;++iSET) { - TrackerHitExtended * trkHitExt = _allSETHits[iSET]; - edm4hep::TrackerHit trkHit = trkHitExt->getTrackerHit(); - int layer = getLayerID(trkHit); - if (layer>=0 && (unsigned)layer < _nLayersSET) - SETHits[layer].push_back(trkHitExt); - } - for (unsigned iL=0; iL< _nLayersSET; ++iL) { // loop over SET layers - TrackerHitExtendedVec hitVec = SETHits[iL]; - int refit = 1; - if(hitVec.empty() == false) AssignOuterHitsToTracks(hitVec,_distCutForSETHits,refit); - } - } + for (int iSET=0;iSET<nSETHits;++iSET) { + TrackerHitExtended * trkHitExt = _allSETHits[iSET]; + edm4hep::TrackerHit trkHit = trkHitExt->getTrackerHit(); + int layer = getLayerID(trkHit); + if (layer>=0 && (unsigned)layer < _nLayersSET) + SETHits[layer].push_back(trkHitExt); + } + for (unsigned iL=0; iL< _nLayersSET; ++iL) { // loop over SET layers + TrackerHitExtendedVec hitVec = SETHits[iL]; + int refit = 1; + if (hitVec.empty() == false) AssignOuterHitsToTracks(hitVec, _distCutForSETHits, refit); + } + } - - CleanUpExtrapolations(); + // Cleaning up helix extrapolations + CleanUpExtrapolations(); + } if (_assignSITHits>0) { // Treatment of left-over SIT hits @@ -5076,7 +5173,18 @@ void FullLDCTrackingAlg::setupGearGeom(){ debug() << " ### gear::SET Parameters from as defined in SSet02 Not Present in GEAR FILE" << endmsg; } } - + + _nLayersETD = 0; + try { + debug() << " filling ETD parameters from gear::ETDParameters " << endmsg; + + const gear::GearParameters& pETDDet = gearMgr->getGearParameters("ETD"); + _nLayersETD = int(pETDDet.getDoubleVals("ETDLayerZ").size()); + } + catch (...) { + debug() << " ### gear::GearParameters ETD Not Present in GEAR FILE" << endmsg; + } + //-- FTD Parameters-- _petalBasedFTDWithOverlaps = false; _nLayersFTD = 0; @@ -5084,7 +5192,7 @@ void FullLDCTrackingAlg::setupGearGeom(){ try{ debug() << " filling FTD parameters from gear::FTDParameters " << endmsg; - const gear::FTDParameters& pFTD =gearMgr->getFTDParameters(); + const gear::FTDParameters& pFTD = gearMgr->getFTDParameters(); const gear::FTDLayerLayout& ftdlayers = pFTD.getFTDLayerLayout() ; _nLayersFTD = ftdlayers.getNLayers() ; diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h index 1bc7d4ae1f35bcca628558fdd5eb3176527e5b75..9864168f820c04a91789be36f09bbd125c2cc496 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h @@ -460,7 +460,7 @@ protected: unsigned int _nLayersSET; - + unsigned int _nLayersETD; // struct FTD_Disk { // int nPetals; // double phi0; @@ -511,10 +511,12 @@ protected: DataHandle<edm4hep::TrackerHitCollection> _FTDPixelTrackerHitColHdl{"FTDPixelTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackerHitCollection> _SITTrackerHitColHdl{"SITSpacePoints", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackerHitCollection> _SETTrackerHitColHdl{"SETSpacePoints", Gaudi::DataHandle::Reader, this}; + DataHandle<edm4hep::TrackerHitCollection> _ETDTrackerHitColHdl{"ETDSpacePoints", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackerHitCollection> _VTXTrackerHitColHdl{"VTXTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackerHitCollection> _inSITRawColHdl{"SITTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackerHitCollection> _inSETRawColHdl{"SETTrackerHits", Gaudi::DataHandle::Reader, this}; + DataHandle<edm4hep::TrackerHitCollection> _inETDRawColHdl{"SETTrackerHits", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackerHitCollection> _inFTDRawColHdl{"FTDStripTrackerHits", Gaudi::DataHandle::Reader, this}; //DataHandle<edm4hep::SimTrackerHitCollection> _inVTXRawColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this}; diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp index e4f4976e187cbf9c3792a5817cccebb193c52bd0..871725fff15b5a5a778432f556c4909c4e442d44 100644 --- a/Service/GearSvc/src/GearSvc.cpp +++ b/Service/GearSvc/src/GearSvc.cpp @@ -712,8 +712,8 @@ StatusCode GearSvc::convertSET(dd4hep::DetElement& set){ n_sensors_per_ladder.push_back(nSensorsPerLadder); setParams->addLayer(nLadders, phi0, supRMin, supOffset, supThickness, supHalfLength, supWidth, 0, senRMin, senOffset, senThickness, senHalfLength, senWidth, 0); } - setParams->setIntVals("n_sensors_per_ladder",n_sensors_per_ladder); - m_gearMgr->setSETParameters( setParams ) ; + setParams->setIntVals("n_sensors_per_ladder", n_sensors_per_ladder); + m_gearMgr->setSETParameters(setParams); return StatusCode::SUCCESS; } diff --git a/Service/TrackSystemSvc/src/MarlinKalTest.cc b/Service/TrackSystemSvc/src/MarlinKalTest.cc index 302675bf183d95f2cd05d7f3b0c793b56f6483b9..bd7fade48e1057c238b0af96b5e643aca9aede7e 100644 --- a/Service/TrackSystemSvc/src/MarlinKalTest.cc +++ b/Service/TrackSystemSvc/src/MarlinKalTest.cc @@ -11,8 +11,10 @@ #include "kaldet/ILDVXDKalDetector.h" #include "kaldet/CEPCVTXKalDetector.h" #include "kaldet/ILDSITKalDetector.h" +#include "kaldet/CEPCITKKalDetector.h" #include "kaldet/ILDSITCylinderKalDetector.h" #include "kaldet/ILDSETKalDetector.h" +#include "kaldet/CEPCOTKKalDetector.h" #include "kaldet/ILDFTDKalDetector.h" #include "kaldet/ILDFTDDiscBasedKalDetector.h" #include "kaldet/ILDTPCKalDetector.h" @@ -127,7 +129,8 @@ namespace MarlinTrk{ bool SIT_found = false ; try{ - ILDSITKalDetector* sitdet = new ILDSITKalDetector( *_gearMgr, _geoSvc ) ; + //ILDSITKalDetector* sitdet = new ILDSITKalDetector( *_gearMgr, _geoSvc ) ; + CEPCITKKalDetector* sitdet = new CEPCITKKalDetector(*_gearMgr, _geoSvc); // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(sitdet); _det->Install( *sitdet ) ; @@ -150,7 +153,8 @@ namespace MarlinTrk{ } try{ - ILDSETKalDetector* setdet = new ILDSETKalDetector( *_gearMgr, _geoSvc ) ; + //ILDSETKalDetector* setdet = new ILDSETKalDetector( *_gearMgr, _geoSvc ); + CEPCOTKKalDetector* setdet = new CEPCOTKKalDetector(*_gearMgr, _geoSvc); // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(setdet); _det->Install( *setdet ) ; diff --git a/Utilities/KalDet/include/kaldet/CEPCITKKalDetector.h b/Utilities/KalDet/include/kaldet/CEPCITKKalDetector.h new file mode 100644 index 0000000000000000000000000000000000000000..7921a442a771564ee6445acf44a0902d633de2bc --- /dev/null +++ b/Utilities/KalDet/include/kaldet/CEPCITKKalDetector.h @@ -0,0 +1,60 @@ +#ifndef __CEPCITKKALDETECTOR +#define __CEPCITKKALDETECTOR__ + +/** Ladder based ITK to be used for ILD DBD studies + * + * @author S.Aplin DESY + */ + +#include "kaltest/TVKalDetector.h" + +#include "TMath.h" + +class TNode; + +namespace gear{ + class GearMgr ; +} + +class IGeomSvc; + +class CEPCITKKalDetector : public TVKalDetector { + +public: + + /** Initialize the ITK from GEAR */ + CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc ); + + +private: + + void setupGearGeom( const gear::GearMgr& gearMgr ) ; + void setupGearGeom( IGeomSvc* geoSvc ); + + int _nLayers ; + double _bZ ; + + bool _isStripDetector; + + struct ITK_Layer { + int nLadders; + int nSensorsPerLadder; + double phi0; + double dphi; + double senRMin; + double supRMin; + double length; + double width; + double offset; + double senThickness; + double supThickness; + double sensorLength; + double stripAngle; + }; + std::vector<ITK_Layer> _ITKgeo; + +}; + + + +#endif diff --git a/Utilities/KalDet/include/kaldet/CEPCOTKKalDetector.h b/Utilities/KalDet/include/kaldet/CEPCOTKKalDetector.h new file mode 100644 index 0000000000000000000000000000000000000000..62aec17960d17548fc0fd4cb5cb984cf220eea14 --- /dev/null +++ b/Utilities/KalDet/include/kaldet/CEPCOTKKalDetector.h @@ -0,0 +1,55 @@ +#ifndef __CEPCOTKKALDETECTOR__ +#define __CEPCOTKKALDETECTOR__ + +/** Ladder based OTK to be used for ILD DBD studies + * + * @author S.Aplin DESY + */ + +#include "kaltest/TVKalDetector.h" + +#include "TMath.h" + +class TNode; + +class IGeomSvc; + +namespace gear{ + class GearMgr ; +} + +class CEPCOTKKalDetector : public TVKalDetector { + public: + /** Initialize the OTK from GEAR */ + CEPCOTKKalDetector(const gear::GearMgr& gearMgr, IGeomSvc* geoSvc=0); + + private: + void setupGearGeom(const gear::GearMgr& gearMgr); + void setupGearGeom(IGeomSvc* geoSvc); + + int _nLayers ; + double _bZ ; + + bool _isStripDetector; + + struct OTK_Layer { + int nLadders; + int nSensorsPerLadder; + double phi0; + double dphi; + double senRMin; + double supRMin; + double senWidth; + double supWidth; + double senOffset; + double supOffset; + double senThickness; + double supThickness; + double senLength; + double supLength; + double sensorLength; + double stripAngle; + }; + std::vector<OTK_Layer> _OTKgeo; +}; +#endif diff --git a/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc b/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc index 30669e42f1823c9ad55892c5ab182b0e078934ab..35c4034da37e7e5eaa483be258c49a25cef2e766 100644 --- a/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc +++ b/Utilities/KalDet/src/ild/common/ILDPlanarMeasLayer.cc @@ -201,12 +201,14 @@ Bool_t ILDPlanarMeasLayer::IsOnSurface(const TVector3 &xx) const #ifdef DEBUG_ISONSURFACE else{ - // streamlog_out(DEBUG4) << "ILDPlanarMeasLayer::IsOnSurface: Point not within boundary x = " << xx.x() << " y = " << xx.y() << " z = " << xx.z() << " r = " << xx.Perp() << " phi = " << xx.Phi() << std::endl; - // streamlog_out(DEBUG4) << "xi = " << xi << " xi_max = " << GetXioffset() + GetXiwidth()/2 << " xi_min = " << GetXioffset() - GetXiwidth()/2 << " zeta = " << zeta << " zeta_min = " << -GetZetawidth()/2 << " zeta_max " << GetZetawidth()/2 << " Xioffset = " << GetXioffset() << std::endl; + std::cout << "ILDPlanarMeasLayer::IsOnSurface: Point not within boundary" + << " x = " << xx.x() << " y = " << xx.y() << " z = " << xx.z() << " r = " << xx.Perp() << " phi = " << xx.Phi() << std::endl; + std::cout << "xi = " << xi << " xi_max = " << GetXioffset() + GetXiwidth()/2 << " xi_min = " << GetXioffset() - GetXiwidth()/2 + << " zeta = " << zeta << " zeta_min = " << -GetZetawidth()/2 << " zeta_max " << GetZetawidth()/2 << " Xioffset = " << GetXioffset() << std::endl; onSurface = false; - // streamlog_out(DEBUG4) << " xi <= GetXioffset() + GetXiwidth()/2 = " << (xi <= GetXioffset() + GetXiwidth()/2) << std::endl ; - // streamlog_out(DEBUG4) << " xi >= GetXioffset() - GetXiwidth()/2 = " << (xi >= GetXioffset() - GetXiwidth()/2) << std::endl; - // streamlog_out(DEBUG4) << " TMath::Abs(zeta) <= GetZetawidth()/2 = " << (TMath::Abs(zeta) <= GetZetawidth()/2) << std::endl; + std::cout << " xi <= GetXioffset() + GetXiwidth()/2 = " << (xi <= GetXioffset() + GetXiwidth()/2) << std::endl; + std::cout << " xi >= GetXioffset() - GetXiwidth()/2 = " << (xi >= GetXioffset() - GetXiwidth()/2) << std::endl; + std::cout << " TMath::Abs(zeta) <= GetZetawidth()/2 = " << (TMath::Abs(zeta) <= GetZetawidth()/2) << std::endl; } #endif @@ -214,28 +216,29 @@ Bool_t ILDPlanarMeasLayer::IsOnSurface(const TVector3 &xx) const #ifdef DEBUG_ISONSURFACE else{ - // streamlog_out(DEBUG4) << "ILDPlanarMeasLayer::IsOnSurface: Point not on surface x = " << xx.x() << " y = " << xx.y() << " z = " << xx.z() << " r = " << xx.Perp() << " phi = " << xx.Phi() << std::endl; - // streamlog_out(DEBUG4) << "Distance from plane " << (xx.X()-GetXc().X())*GetNormal().X() + (xx.Y()-GetXc().Y())*GetNormal().Y() << std::endl; + std::cout << "ILDPlanarMeasLayer::IsOnSurface: Point not on surface" + << " x = " << xx.x() << " y = " << xx.y() << " z = " << xx.z() << " r = " << xx.Perp() << " phi = " << xx.Phi() << std::endl; + std::cout << "Distance from plane " << (xx.X()-GetXc().X())*GetNormal().X() + (xx.Y()-GetXc().Y())*GetNormal().Y() << std::endl; } if( onSurface == false ) { - // streamlog_out(DEBUG) << "x0 " << GetXc().X() << std::endl; - // streamlog_out(DEBUG) << "y0 " << GetXc().Y() << std::endl; - // streamlog_out(DEBUG) << "z0 " << GetXc().Z() << std::endl; - // streamlog_out(DEBUG) << "GetNormal().X() " << GetNormal().X() << std::endl; - // streamlog_out(DEBUG4) << "GetNormal().Y() " << GetNormal().Y() << std::endl; - // streamlog_out(DEBUG4) << "GetNormal().Perp() " << GetNormal().Perp() << std::endl; - // streamlog_out(DEBUG4) << "GetNormal().X()/GetNormal().Perp() " << GetNormal().X()/GetNormal().Perp() << std::endl; - // streamlog_out(DEBUG4) << "GetNormal().Y()/GetNormal().Perp() " << GetNormal().Y()/GetNormal().Perp() << std::endl; - // streamlog_out(DEBUG4) << "xx.X()-GetXc().X() " << xx.X()-GetXc().X() << std::endl; - // streamlog_out(DEBUG4) << "xx.Y()-GetXc().Y() " << xx.Y()-GetXc().Y() << std::endl; - - // streamlog_out(DEBUG4) << "zeta " << zeta << std::endl; - // streamlog_out(DEBUG4) << "xi " << xi << std::endl; - // streamlog_out(DEBUG4) << "zeta half width " << GetZetawidth()/2 << std::endl; - // streamlog_out(DEBUG4) << "xi half width " << GetXiwidth()/2 << std::endl; - // streamlog_out(DEBUG4) << "offset " << GetXioffset() << std::endl; + std::cout << "x0 " << GetXc().X() << std::endl; + std::cout << "y0 " << GetXc().Y() << std::endl; + std::cout << "z0 " << GetXc().Z() << std::endl; + std::cout << "GetNormal().X() " << GetNormal().X() << std::endl; + std::cout << "GetNormal().Y() " << GetNormal().Y() << std::endl; + std::cout << "GetNormal().Perp() " << GetNormal().Perp() << std::endl; + std::cout << "GetNormal().X()/GetNormal().Perp() " << GetNormal().X()/GetNormal().Perp() << std::endl; + std::cout << "GetNormal().Y()/GetNormal().Perp() " << GetNormal().Y()/GetNormal().Perp() << std::endl; + std::cout << "xx.X()-GetXc().X() " << xx.X()-GetXc().X() << std::endl; + std::cout << "xx.Y()-GetXc().Y() " << xx.Y()-GetXc().Y() << std::endl; + + std::cout << "zeta " << zeta << std::endl; + std::cout << "xi " << xi << std::endl; + std::cout << "zeta half width " << GetZetawidth()/2 << std::endl; + std::cout << "xi half width " << GetXiwidth()/2 << std::endl; + std::cout << "offset " << GetXioffset() << std::endl; - // streamlog_out(DEBUG4) << "distance from plane " << (xx.X()-GetXc().X())*GetNormal().X() + (xx.Y()-GetXc().Y())*GetNormal().Y() << std::endl; + std::cout << "distance from plane " << (xx.X()-GetXc().X())*GetNormal().X() + (xx.Y()-GetXc().Y())*GetNormal().Y() << std::endl; } #endif @@ -288,23 +291,25 @@ ILDVTrackHit* ILDPlanarMeasLayer::ConvertLCIOTrkHit(edm4hep::TrackerHit trkhit) dx[1] = trkhit.getCovMatrix(5); bool hit_on_surface = IsOnSurface(hit); - /* - std::cout << "ILDPlanarMeasLayer::ConvertLCIOTrkHit ILDPlanarHit created" - << " for CellID " << trkhit.getCellID() - << " Layer R = " << this->GetXc().Perp() - << " Layer phi = " << this->GetXc().Phi() - << " Layer z0 = " << this->GetXc().Z() - << " u = " << x[0] - << " v = " << x[1] - << " du = " << dx[0] - << " dv = " << dx[1] - << " x = " << hit.x() - << " y = " << hit.y() - << " z = " << hit.z() - << " r = " << hit.Perp() - << " onSurface = " << hit_on_surface - << std::endl ; - */ + +//#define DEBUG_CONVERT 1 +#ifdef DEBUG_CONVERT + std::cout << "ILDPlanarMeasLayer::ConvertLCIOTrkHit ILDPlanarHit created" + << " for CellID " << trkhit.getCellID() + << " Layer R = " << this->GetXc().Perp() + << " Layer phi = " << this->GetXc().Phi() + << " Layer z0 = " << this->GetXc().Z() + << " u = " << x[0] + << " v = " << x[1] + << " du = " << dx[0] + << " dv = " << dx[1] + << " x = " << hit.x() + << " y = " << hit.y() + << " z = " << hit.z() + << " r = " << hit.Perp() + << " onSurface = " << hit_on_surface + << std::endl; +#endif return hit_on_surface ? new ILDPlanarHit( *this , x, dx, this->GetBz(), trkhit) : NULL; } diff --git a/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc b/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc new file mode 100644 index 0000000000000000000000000000000000000000..ffe84e4433f2f9608e04a11f9825b978987b1069 --- /dev/null +++ b/Utilities/KalDet/src/ild/set/CEPCOTKKalDetector.cc @@ -0,0 +1,340 @@ + +#include "kaldet/CEPCOTKKalDetector.h" + +#include "kaldet/MaterialDataBase.h" + +#include "kaldet/ILDParallelPlanarStripMeasLayer.h" + + +#include <UTIL/BitField64.h> +#include <UTIL/ILDConf.h> + +#include "DetInterface/IGeomSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" + +#include <gear/GEAR.h> +#include "gear/BField.h" +#include <gear/ZPlanarParameters.h> +#include <gear/ZPlanarLayerLayout.h> +#include "gearimpl/Util.h" + +#include "CLHEP/Units/SystemOfUnits.h" +#include "DD4hep/DD4hepUnits.h" +#include "TMath.h" + +#include "math.h" +#include <sstream> + +// #include "streamlog/streamlog.h" + +CEPCOTKKalDetector::CEPCOTKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc ) +: TVKalDetector(300) // SJA:FIXME initial size, 300 looks reasonable for ILD, though this would be better stored as a const somewhere +{ + + + // streamlog_out(DEBUG4) << "CEPCOTKKalDetector building OTK detector using GEAR " << std::endl ; + + MaterialDataBase::Instance().registerForService(gearMgr); + + TMaterial & air = *MaterialDataBase::Instance().getMaterial("air"); + TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon"); + TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon"); + + if(geoSvc){ + setupGearGeom( geoSvc ) ; + } + else{ + setupGearGeom( gearMgr ); + } + + if (_isStripDetector) { + // streamlog_out(DEBUG4) << "\t\t building OTK detector as STRIP Detector." << std::endl ; + } else { + // streamlog_out(DEBUG4) << "\t\t building OTK detector as PIXEL Detector." << std::endl ; + } + + + //--The Ladder structure (realistic ladder)-- + int nLadders; + +// Bool_t active = true; + Bool_t dummy = false; + + static const double eps_layer = 1e-6; + static const double eps_sensor = 1e-8; + + UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ; + + for (int layer=0; layer<_nLayers; ++layer) { + nLadders = _OTKgeo[layer].nLadders; + + const double phi0 = _OTKgeo[layer].phi0; + + const double ladder_offset = _OTKgeo[layer].supOffset; + const double ladder_distance = _OTKgeo[layer].supRMin; + const double ladder_thickness = _OTKgeo[layer].supThickness; + const double ladder_width = _OTKgeo[layer].supWidth; + const double ladder_length = _OTKgeo[layer].supLength; + + const double sensitive_offset = _OTKgeo[layer].senOffset; + const double sensitive_distance = _OTKgeo[layer].senRMin; + const double sensitive_thickness = _OTKgeo[layer].senThickness; + const double sensitive_width = _OTKgeo[layer].senWidth; + const double sensitive_length = _OTKgeo[layer].senLength; + + double currPhi; + const double dphi = _OTKgeo[layer].dphi; + + const double stripAngle = pow(-1,layer) *_OTKgeo[layer].stripAngle; + + const int nsensors = _OTKgeo[layer].nSensorsPerLadder; + + const double sensor_length = _OTKgeo[layer].sensorLength; + + for (int ladder=0; ladder<nLadders; ++ladder) { + currPhi = phi0 + (dphi * ladder); + + encoder.reset(); // reset to 0 + + encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::SET; + encoder[lcio::ILDCellID0::side] = 0; + encoder[lcio::ILDCellID0::layer] = layer; + encoder[lcio::ILDCellID0::module] = ladder; + + double z_centre_support = 0.0; + + // check if the sensitive is inside or outside for the support + if (sensitive_distance < ladder_distance) { + + double sen_front_sorting_policy = sensitive_distance + (4 * ladder+0) * eps_layer ; + double sen_back_sorting_policy = sensitive_distance + (4 * ladder+2) * eps_layer ; + double sup_back_sorting_policy = ladder_distance + (4 * ladder+3) * eps_layer ; + + // air - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(air, silicon, sensitive_distance, currPhi, _bZ, sen_front_sorting_policy, sensitive_width, sensitive_length, + sensitive_offset, z_centre_support, sensitive_offset, dummy, -1, "OTKSenFront")); + + for (int isensor=0; isensor<nsensors; ++isensor) { + encoder[lcio::ILDCellID0::sensor] = isensor; + int CellID = encoder.lowWord(); + + double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+1) * eps_layer + eps_sensor * isensor ; + + double z_centre_sensor = -0.5*sensitive_length + (0.5*sensor_length) + (isensor*sensor_length) ; + + if (_isStripDetector) { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarStripMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, + measurement_plane_sorting_policy, sensitive_width, sensor_length, sensitive_offset, + z_centre_sensor, sensitive_offset, stripAngle, CellID, "OTKStripMeaslayer")); + } + else { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, + measurement_plane_sorting_policy, sensitive_width, sensor_length, sensitive_offset, + z_centre_sensor, sensitive_offset, true, CellID, "OTKMeaslayer")); + } + + //std::cout << "CEPCOTKKalDetector add surface with CellID = " << CellID + // << " layer = " << layer << " ladder = " << ladder << " sensor = " << isensor << std::endl; + } + + // sensitive - air boundary + Add(new ILDParallelPlanarMeasLayer(silicon, air, sensitive_distance+sensitive_thickness, currPhi, _bZ, + sen_back_sorting_policy, sensitive_width, sensitive_length, sensitive_offset, + z_centre_support, sensitive_offset, dummy, -1, "OTKSenSupportIntf")); + // air - support boundary, very thin air gap added for different support and sensitive size + Add(new ILDParallelPlanarMeasLayer(air, carbon, ladder_distance+1e6, currPhi, _bZ, sen_back_sorting_policy+eps_layer*0.1, + ladder_width, ladder_length, ladder_offset, z_centre_support, ladder_offset, dummy, -1, "OTKSupSenGap")); + // support - air boundary + Add(new ILDParallelPlanarMeasLayer(carbon, air, ladder_distance+ladder_thickness, currPhi, _bZ, + sup_back_sorting_policy, ladder_width, ladder_length, ladder_offset, + z_centre_support, ladder_offset, dummy,-1,"OTKSupRear")); + } + else { + + double sup_front_sorting_policy = ladder_distance + (4 * ladder+0) * eps_layer ; + double sen_front_sorting_policy = sensitive_distance + (4 * ladder+1) * eps_layer ; + double sen_back_sorting_policy = sensitive_distance + (4 * ladder+3) * eps_layer ; + + // air - support boundary + Add(new ILDParallelPlanarMeasLayer(air, carbon, ladder_distance, currPhi, _bZ, sup_front_sorting_policy, + ladder_width, ladder_length, ladder_offset, z_centre_support, ladder_offset, dummy, -1, "OTKSupFront")); + // support - air boundary, very thin air gap added for different support and sensitive size + Add(new ILDParallelPlanarMeasLayer(carbon, air, sensitive_distance-1e6, currPhi, _bZ, sup_front_sorting_policy+eps_layer*0.1, + ladder_width, ladder_length, ladder_offset, z_centre_support, ladder_offset, dummy, -1, "OTKSupSenGap")); + // air boundary - sensitive + Add(new ILDParallelPlanarMeasLayer(air, silicon, sensitive_distance, currPhi, _bZ, sen_front_sorting_policy, + sensitive_width, sensitive_length, sensitive_offset, z_centre_support, sensitive_offset, dummy, -1, "OTKSenSupportIntf")); + + for (int isensor=0; isensor<nsensors; ++isensor) { + encoder[lcio::ILDCellID0::sensor] = isensor; + int CellID = encoder.lowWord() ; + + double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+2) * eps_layer + eps_sensor * isensor ; + + double z_centre_sensor = -0.5*sensitive_length + (0.5*sensor_length) + (isensor*sensor_length) ; + + if (_isStripDetector) { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarStripMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, + measurement_plane_sorting_policy, sensitive_width, sensor_length, sensitive_offset, + z_centre_sensor, sensitive_offset, stripAngle, CellID, "OTKStripMeaslayer")); + } + else { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, + measurement_plane_sorting_policy, sensitive_width, sensor_length, sensitive_offset, + z_centre_sensor, sensitive_offset, true, CellID, "OTKMeaslayer")); + } + + //std::cout << "CEPCOTKKalDetector add surface with CellID = " << CellID + // << " layer = " << layer << " ladder = " << ladder << " sensor = " << isensor << std::endl ; + } + // sensitive - air boundary + Add(new ILDParallelPlanarMeasLayer(silicon, air, sensitive_distance+sensitive_thickness, currPhi, _bZ, sen_back_sorting_policy, + sensitive_width, sensitive_length, sensitive_offset, z_centre_support, sensitive_offset, dummy, -1, "OTKSenRear")); + } + } + } + + SetOwner(); +} + + + +void CEPCOTKKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ + + const gear::ZPlanarParameters& pOTKDetMain = gearMgr.getSETParameters(); + const gear::ZPlanarLayerLayout& pOTKLayerLayout = pOTKDetMain.getZPlanarLayerLayout(); + + _bZ = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; + + _nLayers = pOTKLayerLayout.getNLayers(); + _OTKgeo.resize(_nLayers); + + bool n_sensors_per_ladder_present = true; + + try { + + std::vector<int> v = pOTKDetMain.getIntVals("n_sensors_per_ladder"); + + } catch (gear::UnknownParameterException& e) { + + n_sensors_per_ladder_present = false; + + } + + double strip_angle_deg = 0.0; + + try { + strip_angle_deg = pOTKDetMain.getDoubleVal("strip_angle_deg"); + } + catch (gear::UnknownParameterException& e) {} + + if (strip_angle_deg==0) _isStripDetector = false; + else _isStripDetector = true; + //SJA:FIXME: for now the support is taken as the same size the sensitive + // if this is not done then the exposed areas of the support would leave a carbon - air boundary, + // which if traversed in the reverse direction to the next boundary then the track would be propagated through carbon + // for a significant distance + //std::cout << "=============OTK strip angle: " << strip_angle_deg << "==============" << std::endl; + for( int layer=0; layer < _nLayers; ++layer){ + _OTKgeo[layer].nLadders = pOTKLayerLayout.getNLadders(layer); + _OTKgeo[layer].phi0 = pOTKLayerLayout.getPhi0(layer); + _OTKgeo[layer].dphi = 2*M_PI / _OTKgeo[layer].nLadders; + _OTKgeo[layer].senRMin = pOTKLayerLayout.getSensitiveDistance(layer); + _OTKgeo[layer].supRMin = pOTKLayerLayout.getLadderDistance(layer); + _OTKgeo[layer].senLength = pOTKLayerLayout.getSensitiveLength(layer)*2.0; // note: gear for historical reasons uses the halflength + _OTKgeo[layer].supLength = pOTKLayerLayout.getLadderLength(layer)*2.0; + _OTKgeo[layer].senWidth = pOTKLayerLayout.getSensitiveWidth(layer); + _OTKgeo[layer].supWidth = pOTKLayerLayout.getLadderWidth(layer); + _OTKgeo[layer].senOffset = pOTKLayerLayout.getSensitiveOffset(layer); + _OTKgeo[layer].supOffset = pOTKLayerLayout.getLadderOffset(layer); + _OTKgeo[layer].senThickness = pOTKLayerLayout.getSensitiveThickness(layer); + _OTKgeo[layer].supThickness = pOTKLayerLayout.getLadderThickness(layer); + + if (n_sensors_per_ladder_present) { + _OTKgeo[layer].nSensorsPerLadder = pOTKDetMain.getIntVals("n_sensors_per_ladder")[layer]; + } + else{ + _OTKgeo[layer].nSensorsPerLadder = 1 ; + } + + _OTKgeo[layer].sensorLength = _OTKgeo[layer].senLength / _OTKgeo[layer].nSensorsPerLadder; + + if (_isStripDetector) { + _OTKgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _OTKgeo[layer].stripAngle = 0.0 ; + } + //std::cout << _OTKgeo[layer].nLadders << " " << _OTKgeo[layer].phi0 << " "<< _OTKgeo[layer].dphi << " " << _OTKgeo[layer].senRMin << " " << _OTKgeo[layer].supRMin << " " + // << _OTKgeo[layer].length << " " << _OTKgeo[layer].width << " " << _OTKgeo[layer].offset << " " << _OTKgeo[layer].senThickness << " " << _OTKgeo[layer].supThickness << " " + // << _OTKgeo[layer].nSensorsPerLadder << " " << _OTKgeo[layer].sensorLength << std::endl; + // streamlog_out(DEBUG0) << " layer = " << layer << std::endl; + // streamlog_out(DEBUG0) << " nSensorsPerLadder = " << _OTKgeo[layer].nSensorsPerLadder << std::endl; + // streamlog_out(DEBUG0) << " sensorLength = " << _OTKgeo[layer].sensorLength << std::endl; + // streamlog_out(DEBUG0) << " stripAngle = " << _OTKgeo[layer].stripAngle << std::endl; + // streamlog_out(DEBUG0) << " _isStripDetector = " << _isStripDetector << std::endl; + + } +} + +void CEPCOTKKalDetector::setupGearGeom( IGeomSvc* geoSvc ){ + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement set; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="OTK") continue; + set = it->second; + } + dd4hep::rec::ZPlanarData* setData = nullptr; + try{ + setData = set.extension<dd4hep::rec::ZPlanarData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << setData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z()/dd4hep::tesla; + + std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& setlayers = setData->layers; + _nLayers = setlayers.size(); + _OTKgeo.resize(_nLayers); + + double strip_angle_deg = setData->angleStrip/CLHEP::degree; + if (strip_angle_deg==0) _isStripDetector = false; + else _isStripDetector = true; + //std::cout << "=============OTK strip angle: " << strip_angle_deg << "==============" << std::endl; + for( int layer=0; layer < _nLayers; ++layer){ + dd4hep::rec::ZPlanarData::LayerLayout& pOTKLayerLayout = setlayers[layer]; + + _OTKgeo[layer].nLadders = pOTKLayerLayout.ladderNumber; + _OTKgeo[layer].phi0 = pOTKLayerLayout.phi0; + _OTKgeo[layer].dphi = 2*M_PI / _OTKgeo[layer].nLadders; + _OTKgeo[layer].senRMin = pOTKLayerLayout.distanceSensitive*CLHEP::cm; + _OTKgeo[layer].supRMin = pOTKLayerLayout.distanceSupport*CLHEP::cm; + _OTKgeo[layer].senLength = pOTKLayerLayout.zHalfSensitive*2.0*CLHEP::cm; // note: gear for historical reasons uses the halflength + _OTKgeo[layer].supLength = pOTKLayerLayout.zHalfSupport*2.0*CLHEP::cm; + _OTKgeo[layer].senWidth = pOTKLayerLayout.widthSensitive*CLHEP::cm; + _OTKgeo[layer].supWidth = pOTKLayerLayout.widthSupport*CLHEP::cm; + _OTKgeo[layer].senOffset = pOTKLayerLayout.offsetSensitive*CLHEP::cm; + _OTKgeo[layer].supOffset = pOTKLayerLayout.offsetSupport*CLHEP::cm; + _OTKgeo[layer].senThickness = pOTKLayerLayout.thicknessSensitive*CLHEP::cm; + _OTKgeo[layer].supThickness = pOTKLayerLayout.thicknessSupport*CLHEP::cm; + _OTKgeo[layer].nSensorsPerLadder = pOTKLayerLayout.sensorsPerLadder; + _OTKgeo[layer].sensorLength = _OTKgeo[layer].senLength / _OTKgeo[layer].nSensorsPerLadder; + + if (_isStripDetector) { + _OTKgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _OTKgeo[layer].stripAngle = 0.0 ; + } + //std::cout << _OTKgeo[layer].nLadders << " " << _OTKgeo[layer].phi0 << " "<< _OTKgeo[layer].dphi << " " << _OTKgeo[layer].senRMin << " " << _OTKgeo[layer].supRMin << " " + // << _OTKgeo[layer].length << " " << _OTKgeo[layer].width << " " << _OTKgeo[layer].offset << " " << _OTKgeo[layer].senThickness << " " << _OTKgeo[layer].supThickness << " " + // << _OTKgeo[layer].nSensorsPerLadder << " " << _OTKgeo[layer].sensorLength << " " << pOTKLayerLayout.lengthSensor << std::endl; + } +} diff --git a/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc b/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc new file mode 100644 index 0000000000000000000000000000000000000000..ae897cc8adb9368be0c8e8b1658ac1528610aa5d --- /dev/null +++ b/Utilities/KalDet/src/ild/sit/CEPCITKKalDetector.cc @@ -0,0 +1,320 @@ + +#include "kaldet/CEPCITKKalDetector.h" + +#include "MaterialDataBase.h" + +#include "ILDParallelPlanarStripMeasLayer.h" + + +#include <UTIL/BitField64.h> +#include <UTIL/ILDConf.h> + +#include "DetInterface/IGeomSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Units/SystemOfUnits.h" +#include "DD4hep/DD4hepUnits.h" + +#include <gear/GEAR.h> +#include "gear/BField.h" +#include <gear/ZPlanarParameters.h> +#include <gear/ZPlanarLayerLayout.h> +#include "gearimpl/Util.h" + +#include "TMath.h" + +#include "math.h" +#include <sstream> + +// #include "streamlog/streamlog.h" + +CEPCITKKalDetector::CEPCITKKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc ) +: TVKalDetector(300) // SJA:FIXME initial size, 300 looks reasonable for ILD, though this would be better stored as a const somewhere +{ + // std::cout << "CEPCITKKalDetector building ITK detector using GEAR " << std::endl ; + + MaterialDataBase::Instance().registerForService(gearMgr, geoSvc); + + TMaterial & air = *MaterialDataBase::Instance().getMaterial("air"); + TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon"); + TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("carbon"); + + if(geoSvc){ + this->setupGearGeom(geoSvc); + } + else{ + this->setupGearGeom(gearMgr) ; + } + + if (_isStripDetector) { + // streamlog_out(DEBUG4) << "\t\t building ITK detector as STRIP Detector." << std::endl ; + } else { + // streamlog_out(DEBUG4) << "\t\t building ITK detector as PIXEL Detector." << std::endl ; + } + + //--The Ladder structure (realistic ladder)-- + int nLadders; + + // Bool_t active = true; + Bool_t dummy = false; + + static const double eps_layer = 1e-6; + static const double eps_sensor = 1e-8; + + UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ; + + for (int layer=0; layer<_nLayers; ++layer) { + + double offset = _ITKgeo[layer].offset ; + double xioffset = 0.0; + double z_centre_support = 0.0; + + nLadders = _ITKgeo[layer].nLadders ; + + const double phi0 = _ITKgeo[layer].phi0 ; + + const double ladder_distance = _ITKgeo[layer].supRMin ; + const double ladder_thickness = _ITKgeo[layer].supThickness ; + + const double sensitive_distance = _ITKgeo[layer].senRMin ; + const double sensitive_thickness = _ITKgeo[layer].senThickness ; + + const double width = _ITKgeo[layer].width ; + const double length = _ITKgeo[layer].length; + + double currPhi; + const double dphi = _ITKgeo[layer].dphi ; + + const double stripAngle = pow(-1,layer) *_ITKgeo[layer].stripAngle; + + const int nsensors = _ITKgeo[layer].nSensorsPerLadder; + + const double sensor_length = _ITKgeo[layer].sensorLength; + + // TODO: sorting_policy for overlap region like VXD + for (int ladder=0; ladder<nLadders; ++ladder) { + + currPhi = phi0 + (dphi * ladder); + + encoder.reset() ; // reset to 0 + + encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::SIT; + encoder[lcio::ILDCellID0::side] = 0 ; + encoder[lcio::ILDCellID0::layer] = layer ; + encoder[lcio::ILDCellID0::module] = ladder ; + + // check if the sensitive is inside or outside for the support + if( sensitive_distance < ladder_distance ) { + + double sen_front_sorting_policy = sensitive_distance + (4 * ladder+0) * eps_layer ; + double sen_back_sorting_policy = sensitive_distance + (4 * ladder+2) * eps_layer ; + double sup_back_sorting_policy = ladder_distance + (4 * ladder+3) * eps_layer ; + + // air - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(air, silicon, sensitive_distance, currPhi, _bZ, sen_front_sorting_policy, width, length, offset, z_centre_support, offset, dummy,-1,"ITKSenFront")) ; + + for (int isensor=0; isensor<nsensors; ++isensor) { + + encoder[lcio::ILDCellID0::sensor] = isensor ; + int CellID = encoder.lowWord() ; + + double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+1) * eps_layer + eps_sensor * isensor ; + + double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor*sensor_length) ; + + if (_isStripDetector) { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarStripMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, measurement_plane_sorting_policy, width, sensor_length, offset, z_centre_sensor, offset, stripAngle, CellID, "ITKStripMeaslayer")) ; + } + else { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, measurement_plane_sorting_policy, width, sensor_length, offset, z_centre_sensor, offset, true, CellID, "ITKMeaslayer")) ; + } + + //std::cout << "CEPCITKKalDetector add surface with CellID = " << CellID << std::endl ; + } + + // sensitive - support boundary + Add(new ILDParallelPlanarMeasLayer(silicon, carbon, sensitive_distance+sensitive_thickness, currPhi, _bZ, sen_back_sorting_policy, width, length, offset, z_centre_support, offset, dummy,-1,"ITKSenSupportIntf" )) ; + + // support - air boundary + Add(new ILDParallelPlanarMeasLayer(carbon, air, ladder_distance+ladder_thickness, currPhi, _bZ, sup_back_sorting_policy, width, length, offset, z_centre_support, offset, dummy,-1,"ITKSupRear" )) ; + } + else { + + double sup_front_sorting_policy = ladder_distance + (4 * ladder+0) * eps_layer ; + double sen_front_sorting_policy = sensitive_distance + (4 * ladder+1) * eps_layer ; + double sen_back_sorting_policy = sensitive_distance + (4 * ladder+3) * eps_layer ; + + // air - support boundary + Add(new ILDParallelPlanarMeasLayer(air, carbon, ladder_distance, currPhi, _bZ, sup_front_sorting_policy, width, length, offset, z_centre_support, offset, dummy,-1,"ITKSupFront")) ; + + // support boundary - sensitive + Add(new ILDParallelPlanarMeasLayer(carbon, silicon, sensitive_distance, currPhi, _bZ, sen_front_sorting_policy, width, length, offset, z_centre_support, offset, dummy,-1,"ITKSenSupportIntf" )) ; + + for (int isensor=0; isensor<nsensors; ++isensor) { + + encoder[lcio::ILDCellID0::sensor] = isensor ; + int CellID = encoder.lowWord() ; + + double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+2) * eps_layer + eps_sensor * isensor ; + + double z_centre_sensor = -0.5*length + (0.5*sensor_length) + (isensor*sensor_length) ; + + if (_isStripDetector) { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarStripMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, measurement_plane_sorting_policy, width, sensor_length, offset, z_centre_sensor, offset, stripAngle, CellID, "ITKStripMeaslayer")) ; + } else { + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*0.5, currPhi, _bZ, measurement_plane_sorting_policy, width, sensor_length, offset, z_centre_sensor, offset, true, CellID, "ITKMeaslayer")) ; + } + + //std::cout << "CEPCITKKalDetector add surface with CellID = " << CellID << std::endl ; + } + + // support - air boundary + Add(new ILDParallelPlanarMeasLayer(silicon, air, sensitive_distance+sensitive_thickness, currPhi, _bZ, sen_back_sorting_policy, width, length, offset, z_centre_support, offset, dummy,-1,"ITKSenRear" )) ; + } + } + } + + SetOwner(); +} + + + +void CEPCITKKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ + + const gear::ZPlanarParameters& pITKDetMain = gearMgr.getSITParameters(); + const gear::ZPlanarLayerLayout& pITKLayerLayout = pITKDetMain.getZPlanarLayerLayout(); + + _bZ = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; + + _nLayers = pITKLayerLayout.getNLayers(); + _ITKgeo.resize(_nLayers); + + bool n_sensors_per_ladder_present = true; + + try { + + std::vector<int> v = pITKDetMain.getIntVals("n_sensors_per_ladder"); + + } catch (gear::UnknownParameterException& e) { + + n_sensors_per_ladder_present = false; + + } + + double strip_angle_deg = 0.0; + _isStripDetector = false; + + try { + + strip_angle_deg = pITKDetMain.getDoubleVal("strip_angle_deg"); + if (strip_angle_deg !=0 ) _isStripDetector = true; + + } catch (gear::UnknownParameterException& e) { + + _isStripDetector = false; + + } + + + //SJA:FIXME: for now the support is taken as the same size the sensitive + // if this is not done then the exposed areas of the support would leave a carbon - air boundary, + // which if traversed in the reverse direction to the next boundary then the track would be propagated through carbon + // for a significant distance + //std::cout << "=============ITK strip angle: " << strip_angle_deg << "==============" << std::endl; + for( int layer=0; layer < _nLayers; ++layer){ + + _ITKgeo[layer].nLadders = pITKLayerLayout.getNLadders(layer); + _ITKgeo[layer].phi0 = pITKLayerLayout.getPhi0(layer); + _ITKgeo[layer].dphi = 2*M_PI / _ITKgeo[layer].nLadders; + _ITKgeo[layer].senRMin = pITKLayerLayout.getSensitiveDistance(layer); + _ITKgeo[layer].supRMin = pITKLayerLayout.getLadderDistance(layer); + _ITKgeo[layer].length = pITKLayerLayout.getSensitiveLength(layer)*2.0; // note: gear for historical reasons uses the halflength + _ITKgeo[layer].width = pITKLayerLayout.getSensitiveWidth(layer); + _ITKgeo[layer].offset = pITKLayerLayout.getSensitiveOffset(layer); + _ITKgeo[layer].senThickness = pITKLayerLayout.getSensitiveThickness(layer); + _ITKgeo[layer].supThickness = pITKLayerLayout.getLadderThickness(layer); + + if (n_sensors_per_ladder_present) { + _ITKgeo[layer].nSensorsPerLadder = pITKDetMain.getIntVals("n_sensors_per_ladder")[layer]; + } + else{ + _ITKgeo[layer].nSensorsPerLadder = 1 ; + } + + _ITKgeo[layer].sensorLength = _ITKgeo[layer].length / _ITKgeo[layer].nSensorsPerLadder; + + + if (_isStripDetector) { + _ITKgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _ITKgeo[layer].stripAngle = 0.0 ; + } + //std::cout << _ITKgeo[layer].nLadders << " " << _ITKgeo[layer].phi0 << " "<< _ITKgeo[layer].dphi << " " + // << _ITKgeo[layer].senRMin << " " << _ITKgeo[layer].supRMin << " " << _ITKgeo[layer].length << " " + // << _ITKgeo[layer].width << " " << _ITKgeo[layer].offset << " " << _ITKgeo[layer].senThickness << " " + // << _ITKgeo[layer].supThickness << " " << _ITKgeo[layer].nSensorsPerLadder << " " << _ITKgeo[layer].sensorLength << std::endl; + } + //std::cout << " _isStripDetector = " << _isStripDetector << std::endl; +} + +void CEPCITKKalDetector::setupGearGeom( IGeomSvc* geoSvc ){ + + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement sit; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="ITK") continue; + sit = it->second; + } + dd4hep::rec::ZPlanarData* sitData = nullptr; + try{ + sitData = sit.extension<dd4hep::rec::ZPlanarData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << sitData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z()/dd4hep::tesla; + + std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& sitlayers = sitData->layers; + _nLayers = sitlayers.size(); + _ITKgeo.resize(_nLayers); + + _isStripDetector = false; + double strip_angle_deg = sitData->angleStrip/CLHEP::degree; + if(strip_angle_deg!=0){ + _isStripDetector = true; + } + //std::cout << "=============ITK strip angle: " << strip_angle_deg << "==============" << std::endl; + for( int layer=0; layer < _nLayers; ++layer){ + dd4hep::rec::ZPlanarData::LayerLayout& pITKLayerLayout = sitlayers[layer]; + + _ITKgeo[layer].nLadders = pITKLayerLayout.ladderNumber; + _ITKgeo[layer].phi0 = pITKLayerLayout.phi0; + _ITKgeo[layer].dphi = 2*M_PI / _ITKgeo[layer].nLadders; + _ITKgeo[layer].senRMin = pITKLayerLayout.distanceSensitive*CLHEP::cm; + _ITKgeo[layer].supRMin = pITKLayerLayout.distanceSupport*CLHEP::cm; + _ITKgeo[layer].length = pITKLayerLayout.zHalfSensitive*2.0*CLHEP::cm; // note: gear for historical reasons uses the halflength + _ITKgeo[layer].width = pITKLayerLayout.widthSensitive*CLHEP::cm; + _ITKgeo[layer].offset = pITKLayerLayout.offsetSensitive*CLHEP::cm; + _ITKgeo[layer].senThickness = pITKLayerLayout.thicknessSensitive*CLHEP::cm; + _ITKgeo[layer].supThickness = pITKLayerLayout.thicknessSupport*CLHEP::cm; + _ITKgeo[layer].nSensorsPerLadder = pITKLayerLayout.sensorsPerLadder; + _ITKgeo[layer].sensorLength = _ITKgeo[layer].length / _ITKgeo[layer].nSensorsPerLadder; + + if (_isStripDetector) { + _ITKgeo[layer].stripAngle = strip_angle_deg * M_PI/180 ; + } else { + _ITKgeo[layer].stripAngle = 0.0 ; + } + //std::cout << _ITKgeo[layer].nLadders << " " << _ITKgeo[layer].phi0 << " "<< _ITKgeo[layer].dphi << " " << _ITKgeo[layer].senRMin << " " << _ITKgeo[layer].supRMin << " " + // << _ITKgeo[layer].length << " " << _ITKgeo[layer].width << " " << _ITKgeo[layer].offset << " " << _ITKgeo[layer].senThickness << " " << _ITKgeo[layer].supThickness << " " + // << _ITKgeo[layer].nSensorsPerLadder << " " << _ITKgeo[layer].sensorLength << " " << pITKLayerLayout.lengthSensor << std::endl; + } +} diff --git a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc index 3b7528b5258cf58a02fb428c4839c9f59e057b5f..093b6bdabb6a058dd2ce5ac6c1024a275805dd1c 100644 --- a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc +++ b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc @@ -1,5 +1,5 @@ -#include "CEPCVTXKalDetector.h" +#include "kaldet/CEPCVTXKalDetector.h" #include "MaterialDataBase.h" diff --git a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.h b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.h deleted file mode 100644 index 7184f78c5c8b25ab1b2373fa3be15cf7f36a5b38..0000000000000000000000000000000000000000 --- a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.h +++ /dev/null @@ -1,65 +0,0 @@ -#ifndef __CEPCVTXKALDETECTOR__ -#define __CEPCVTXKALDETECTOR__ - -/** Ladder based VXD to be used for ILD DBD studies - * - * @author S.Aplin DESY - */ - -#include "kaltest/TVKalDetector.h" - -#include "TMath.h" - -class TNode; -class IGeomSvc; - -namespace gear{ - class GearMgr ; -} - - -class CEPCVTXKalDetector : public TVKalDetector { - -public: - - CEPCVTXKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc); - - -private: - - void setupGearGeom( const gear::GearMgr& gearMgr ); - void setupGearGeom( IGeomSvc* geoSvc) ; - - int _nLayers[2]; - double _bZ; - - double _relative_position_of_measurement_surface; - - struct VXD_Layer { - int nLadders; - double phi0; - double dphi; - double senRMin; - double supRMin; - double length; - double width; - double offset; - double senThickness; - double supThickness; - }; - std::vector<VXD_Layer> _VXDgeo; - - struct STT_Layer { - int id; - double length; - double senRMin; - double senThickness; - double supRMin; - double supThickness; - double phi0; - double rgap; - double dphi; - }; - std::vector<STT_Layer> _STTgeo; -}; -#endif