diff --git a/CMakeLists.txt b/CMakeLists.txt index aedbe21b4..08c8bf3d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,7 @@ ENDIF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) option(BUILD_TESTING "Enable and build tests" ON) option(INSTALL_COMPACT_FILES "Copy compact files to install area" OFF) +option(INSTALL_BEAMPIPE_STL_FILES "Download CAD files for building the detailed beampipe" OFF) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -68,6 +69,7 @@ endif() file(GLOB sources ./detector/tracker/*.cpp ./detector/calorimeter/*.cpp + ./detector/calorimeter/dual-readout/src/*.cpp ./detector/fcal/*.cpp ./detector/muonSystem/*.cpp ./detector/other/*.cpp @@ -107,9 +109,15 @@ if(NOT DCH_INFO_H_EXIST) message(WARNING "Subdetector ${FILES_DEPENDINGON_DCH_INFO_H} will not be built because header file DDRec/DCH_info.h was not found") endif() +find_package(EDM4HEP) file(GLOB G4sources ./plugins/TPCSDAction.cpp ./plugins/CaloPreShowerSDAction.cpp + ./plugins/FiberDRCaloSDAction.h + ./plugins/FiberDRCaloSDAction.cpp + ./plugins/Geant4Output2EDM4hep_DRC.cpp + ./plugins/DRCaloFastSimModel.cpp + ./plugins/DRCaloFastSimModel.h ) if(DD4HEP_USE_PYROOT) @@ -125,8 +133,11 @@ add_library(lcgeo ALIAS k4geo) target_include_directories(${PackageName} PRIVATE ${PROJECT_SOURCE_DIR}/detector/include ) target_include_directories(${PackageName}G4 PRIVATE ${PROJECT_SOURCE_DIR}/detector/include ) +target_include_directories(${PackageName} PRIVATE ${PROJECT_SOURCE_DIR}/detector/calorimeter/dual-readout/include ) +target_include_directories(${PackageName}G4 PRIVATE ${PROJECT_SOURCE_DIR}/detector/calorimeter/dual-readout/include ) + target_link_libraries(${PackageName} DD4hep::DDCore DD4hep::DDRec DD4hep::DDParsers ROOT::Core detectorSegmentations) -target_link_libraries(${PackageName}G4 DD4hep::DDCore DD4hep::DDRec DD4hep::DDParsers DD4hep::DDG4 ROOT::Core ${Geant4_LIBRARIES}) +target_link_libraries(${PackageName}G4 DD4hep::DDCore DD4hep::DDRec DD4hep::DDParsers DD4hep::DDG4 ROOT::Core EDM4HEP::edm4hep EDM4HEP::edm4hepDict podio::podio podio::podioDict podio::podioRootIO ${Geant4_LIBRARIES}) if(K4GEO_USE_LCIO) target_link_libraries(${PackageName} LCIO::lcio) @@ -184,3 +195,47 @@ write_basic_package_version_file( install(FILES ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_CMAKEDIR}/k4geoConfig.cmake ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_CMAKEDIR}/k4geoConfigVersion.cmake DESTINATION ${CMAKE_INSTALL_CMAKEDIR} ) + +if(INSTALL_BEAMPIPE_STL_FILES) + + set(STL_FILES + "AlBeMet162_30042024.stl" + "Copper_pipe_28092023.stl" + "Gold_19042024.stl" + "Paraffine_19042024.stl" + "Tungsten_mask_02102023.stl" + "Water_30042024.stl" + ) + # Set main FCC url + set(FCC_URL "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco") + set(STL_PATH "MDI_o1_v01/stl_files/Pipe_240430") + + # Set the output directory where the file will be placed + set(OUTPUT_DIR "${PROJECT_SOURCE_DIR}/FCCee/MDI/compact/${STL_PATH}") + file(MAKE_DIRECTORY ${OUTPUT_DIR}) + + foreach(STL_FILE ${STL_FILES}) + set(FULL_URL "${FCC_URL}/MDI/${STL_PATH}/${STL_FILE}") + message(DEBUG "Downloading file ${FULL_URL}") + set(OUTPUT_FILE "${OUTPUT_DIR}/${STL_FILE}") + + if(EXISTS "${OUTPUT_FILE}") + message(STATUS "File ${STL_FILE} already exists. Skipping download.") + else() + # Download the file + file(DOWNLOAD ${FULL_URL} ${OUTPUT_FILE} + SHOW_PROGRESS + STATUS download_status) + + list(GET download_status 0 status_code) + if(NOT status_code EQUAL 0) + list(GET download_status 1 error_message) + message(FATAL_ERROR "Failed to download file: ${error_message}") + endif() + endif() + endforeach() + + file(MAKE_DIRECTORY share/k4geo/FCCee) + INSTALL(DIRECTORY ${OUTPUT_DIR} DESTINATION share/k4geo/FCCee/MDI/compact/MDI_o1_v01/stl_files ) + +endif() diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml index 4858ee905..7606b7f23 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml @@ -34,6 +34,8 @@ + @@ -42,9 +44,9 @@ - - + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/DectDimensions.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/DectDimensions.xml index c09dbea21..d100e0a0f 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/DectDimensions.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/DectDimensions.xml @@ -103,6 +103,7 @@ + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged.xml index 7a978c8d9..5b7dc11fa 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged.xml @@ -88,13 +88,13 @@ - + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 - + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_calibration.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_calibration.xml index e3e45f658..c8c79f253 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_calibration.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_calibration.xml @@ -80,13 +80,13 @@ - + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 - + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_upstream.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_upstream.xml index 31ea25fa0..b4e05e073 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_upstream.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalBarrel_thetamodulemerged_upstream.xml @@ -80,13 +80,13 @@ - + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 - + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 @@ -154,4 +154,4 @@ - \ No newline at end of file + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml new file mode 100644 index 000000000..547053962 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml @@ -0,0 +1,133 @@ + + + + + + Liquid argon EM calorimeter endcap design. + Electromagnetic part (EMEC) includes lead+steel absorber. + Turbine geometry. + + + + + + + + + + + + + + + + + + + + + + offset + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,cryo:1,type:3,subtype:3,side:-2,wheel:3,layer:8,module:17,rho:8,z:8 + + + --> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml new file mode 100644 index 000000000..3d07d3671 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml @@ -0,0 +1,125 @@ + + + + + Liquid argon EM calorimeter endcap design. + Electromagnetic part (EMEC) includes lead+steel absorber. + Turbine geometry. + + + + + + + + + + + + + + + + + + + + + + offset + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,cryo:1,type:3,subtype:3,side:-2,wheel:3,layer:8,module:17,rho:8,z:8 + + + --> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/FCCee/ALLEGRO/compact/README.md b/FCCee/ALLEGRO/compact/README.md index 2d8dd9f66..6225008c8 100644 --- a/FCCee/ALLEGRO/compact/README.md +++ b/FCCee/ALLEGRO/compact/README.md @@ -8,4 +8,5 @@ Known caveat: the drift chamber has a larger z extent than in the IDEA detector ALLEGRO_o1_v03: with respect to v02 it features an ECal barrel with 11 layers and cell corners projective along phi. The vertex detector and drift chamber are now taken directly from IDEA_o1_v03, this effectively updates both the vertex detector (which was taken from an old CLD version) and the drift chamber (which was corresponding to IDEA_o1_v02/DriftChamber_o1_v01.xml). The z-extent of the drift chamber is now unchanged w.r.t. the IDEA detector (2 m) since it requires optimization anyway. -Magnetic fields (solenoid + MDI) have been added. \ No newline at end of file +Magnetic fields (solenoid + MDI) have been added. +Added "turbine-style" endcap ecal, and invoke this in the top-level xml (replacing the coneCyro geometry). diff --git a/FCCee/CLD/compact/CLD_o2_v06/CLD_o2_v06.xml b/FCCee/CLD/compact/CLD_o2_v06/CLD_o2_v06.xml index 853cd02d4..4599dc34d 100644 --- a/FCCee/CLD/compact/CLD_o2_v06/CLD_o2_v06.xml +++ b/FCCee/CLD/compact/CLD_o2_v06/CLD_o2_v06.xml @@ -2,12 +2,12 @@ xmlns:xs="http://www.w3.org/2001/XMLSchema" xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> - + version="6"> The compact format for the FCCee Detector design diff --git a/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml b/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml index cf080fa0c..e662f7a58 100644 --- a/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml +++ b/FCCee/IDEA/compact/IDEA_o1_v03/DectDimensions_IDEA_o1_v03.xml @@ -53,7 +53,10 @@ - + + + + @@ -107,45 +110,6 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -217,6 +181,20 @@ + + + + + + + + + + + + + + @@ -231,6 +209,7 @@ + @@ -247,6 +226,8 @@ + + @@ -279,6 +260,18 @@ + + + + + + + + + + + + diff --git a/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml b/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml new file mode 100644 index 000000000..8961ffac4 --- /dev/null +++ b/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml @@ -0,0 +1,698 @@ + + + + + + + + The compact format for the dual-readout calorimeter (for FCCee IDEA) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,assembly:1,eta:-8,phi:9,x:32:-11,y:-9,c:1,module:2 + + + + diff --git a/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml b/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml index 66a5c49bd..949c55eb7 100644 --- a/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml +++ b/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml @@ -37,7 +37,10 @@ + + @@ -57,6 +60,10 @@ + + + + diff --git a/FCCee/IDEA/compact/README.md b/FCCee/IDEA/compact/README.md index 192d63f55..fed6cf941 100644 --- a/FCCee/IDEA/compact/README.md +++ b/FCCee/IDEA/compact/README.md @@ -14,4 +14,9 @@ Based on o1_v01 but with a detailed description of the vertex detector, drift ch IDEA_o1_v03 ------------ -Based on o1_v02 but replacing the drift chamber (o1, v01) for the lightweight implementation based on twisted tubes (o1, v02). NB: production threshold and step limit physics have to be tuned for the drift chamber. July 2024: Added a detailed version of the muon system. +Based on o1_v02 but replacing the drift chamber (o1, v01) for the lightweight implementation based on twisted tubes (o1, +v02). NB: production threshold and step limit physics have to be tuned for the drift chamber. + +July 2024: Added a detailed version of the muon system. + +The monolithic fiber dual-readout calorimeter (o1, v01) is added to the directory, but commented in the main file IDEA_o1_v03.xml for the sake of speed. Please remove the comments to include this calorimeter in the full IDEA detector. diff --git a/FCCee/MDI/compact/MDI_o1_v01/Beampipe_CADimport_o1_v02.xml b/FCCee/MDI/compact/MDI_o1_v01/Beampipe_CADimport_o1_v02.xml index c4a3132f6..21df662c8 100644 --- a/FCCee/MDI/compact/MDI_o1_v01/Beampipe_CADimport_o1_v02.xml +++ b/FCCee/MDI/compact/MDI_o1_v01/Beampipe_CADimport_o1_v02.xml @@ -19,7 +19,7 @@ - @@ -31,7 +31,7 @@ - @@ -44,7 +44,7 @@ - @@ -57,7 +57,7 @@ - @@ -69,7 +69,7 @@ - @@ -81,7 +81,7 @@ - diff --git a/FCCee/MDI/compact/README.md b/FCCee/MDI/compact/README.md index 0b177c6d5..13222dae2 100644 --- a/FCCee/MDI/compact/README.md +++ b/FCCee/MDI/compact/README.md @@ -5,9 +5,10 @@ aciarma - 08/07/24 -- BeamInstrumentation_o1_v01.xml : compensating and screening solenoids -- HOMAbsorber.xml : old model of the HOM absorbers. Not needed anymore with new low impedance beam pipe. -- MDI_o1_v01 (STL FILES NOT IN THIS REPO!!! WILL BE ADDED LATER) +- MDI_o1_v01 -- Beampipe_CADimport_o1_v02.xml : import CAD models for engineered beam pipe (by F. Fransesini/INFN-LNF) -These .stl files are hosted [here](https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/MDI/MDI_o1_v01/) +These .stl files are hosted [here](https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/MDI/MDI_o1_v01/). +The CMake option `INSTALL_BEAMPIPE_STL_FILES=ON` downloads these STL. -- stl_files/Pipe_240430 ├── AlBeMet162_30042024.stl : central and elliptoconical chambers, with cooling manifolds ├── Copper_pipe_28092023.stl : low impedance beam pipe separation region @@ -15,4 +16,4 @@ These .stl files are hosted [here](https://fccsw.web.cern.ch/fccsw/filesForSimDi ├── Paraffine_19042024.stl : cooling for central chamber ├── Tungsten_mask_02102023.stl : SR masks 2.1m upstream └── Water_30042024.stl : cooling for elliptoconical chambers --- BeamInstrumentation_o1_v01.xml : compensating and screening solenoids \ No newline at end of file +-- BeamInstrumentation_o1_v01.xml : compensating and screening solenoids diff --git a/ILD/compact/ILD_l5_v11/InnerTrackerILD_o1_v01_00.xml b/ILD/compact/ILD_l5_v11/InnerTrackerILD_o1_v01_00.xml index 1d5142850..8f1df743f 100644 --- a/ILD/compact/ILD_l5_v11/InnerTrackerILD_o1_v01_00.xml +++ b/ILD/compact/ILD_l5_v11/InnerTrackerILD_o1_v01_00.xml @@ -120,13 +120,8 @@ adapted by DJeans to fit into ILD model at FCCee - - diff --git a/detector/CaloTB/CaloPrototype_v01.cpp b/detector/CaloTB/CaloPrototype_v01.cpp index 47fc83585..18816b815 100644 --- a/detector/CaloTB/CaloPrototype_v01.cpp +++ b/detector/CaloTB/CaloPrototype_v01.cpp @@ -8,11 +8,7 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Layering.h" #include "XML/Utilities.h" -#include "DDRec/DetectorData.h" -#include "LcgeoExceptions.h" -#include -#include using namespace std; diff --git a/detector/CaloTB/CaloPrototype_v02.cpp b/detector/CaloTB/CaloPrototype_v02.cpp index d443f2803..d1bc40ee0 100644 --- a/detector/CaloTB/CaloPrototype_v02.cpp +++ b/detector/CaloTB/CaloPrototype_v02.cpp @@ -8,12 +8,8 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Layering.h" #include "XML/Utilities.h" -#include "DDRec/DetectorData.h" #include "DDSegmentation/TiledLayerGridXY.h" -#include "LcgeoExceptions.h" -#include -#include using namespace std; diff --git a/detector/CaloTB/CaloPrototype_v03.cpp b/detector/CaloTB/CaloPrototype_v03.cpp index 6484accc1..360487aea 100644 --- a/detector/CaloTB/CaloPrototype_v03.cpp +++ b/detector/CaloTB/CaloPrototype_v03.cpp @@ -10,7 +10,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Layering.h" #include "XML/Utilities.h" -#include "LcgeoExceptions.h" #include diff --git a/detector/CaloTB/TBecal4d.cpp b/detector/CaloTB/TBecal4d.cpp index 5765ec086..168aa00ac 100644 --- a/detector/CaloTB/TBecal4d.cpp +++ b/detector/CaloTB/TBecal4d.cpp @@ -8,11 +8,7 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Layering.h" #include "XML/Utilities.h" -#include "DDRec/DetectorData.h" -#include "DDSegmentation/TiledLayerGridXY.h" -#include "LcgeoExceptions.h" -#include #include using namespace std; diff --git a/detector/CaloTB/TBhcal4d.cpp b/detector/CaloTB/TBhcal4d.cpp index 2d2e8c6cb..39c462717 100644 --- a/detector/CaloTB/TBhcal4d.cpp +++ b/detector/CaloTB/TBhcal4d.cpp @@ -8,11 +8,7 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Layering.h" #include "XML/Utilities.h" -#include "DDRec/DetectorData.h" -#include "DDSegmentation/TiledLayerGridXY.h" -#include "LcgeoExceptions.h" -#include #include using namespace std; diff --git a/detector/calorimeter/ECalBarrel_o1_v01_geo.cpp b/detector/calorimeter/ECalBarrel_o1_v01_geo.cpp index f8b8266a7..eb9da2a67 100644 --- a/detector/calorimeter/ECalBarrel_o1_v01_geo.cpp +++ b/detector/calorimeter/ECalBarrel_o1_v01_geo.cpp @@ -12,9 +12,7 @@ //========================================================================= #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/Segmentation.h" diff --git a/detector/calorimeter/ECalBarrel_o1_v02_geo.cpp b/detector/calorimeter/ECalBarrel_o1_v02_geo.cpp index 1274bb77d..6fc9f3422 100644 --- a/detector/calorimeter/ECalBarrel_o1_v02_geo.cpp +++ b/detector/calorimeter/ECalBarrel_o1_v02_geo.cpp @@ -15,9 +15,7 @@ #define VERBOSE_LEVEL 2 #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/Segmentation.h" diff --git a/detector/calorimeter/ECalBarrel_o1_v03_geo.cpp b/detector/calorimeter/ECalBarrel_o1_v03_geo.cpp index aaa8540ca..97b154e6f 100644 --- a/detector/calorimeter/ECalBarrel_o1_v03_geo.cpp +++ b/detector/calorimeter/ECalBarrel_o1_v03_geo.cpp @@ -16,9 +16,7 @@ #define VERBOSE_LEVEL 0 #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/Segmentation.h" diff --git a/detector/calorimeter/ECalBarrel_o2_v03_geo.cpp b/detector/calorimeter/ECalBarrel_o2_v03_geo.cpp index 20e24f26d..14e12c102 100644 --- a/detector/calorimeter/ECalBarrel_o2_v03_geo.cpp +++ b/detector/calorimeter/ECalBarrel_o2_v03_geo.cpp @@ -16,9 +16,7 @@ #define VERBOSE_LEVEL 0 #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/Segmentation.h" diff --git a/detector/calorimeter/ECalBarrel_o2_v04_geo.cpp b/detector/calorimeter/ECalBarrel_o2_v04_geo.cpp index 40d57f439..024cd9c9a 100644 --- a/detector/calorimeter/ECalBarrel_o2_v04_geo.cpp +++ b/detector/calorimeter/ECalBarrel_o2_v04_geo.cpp @@ -17,9 +17,7 @@ #define VERBOSE_LEVEL 0 #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/Segmentation.h" diff --git a/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp b/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp new file mode 100644 index 000000000..51d79a656 --- /dev/null +++ b/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp @@ -0,0 +1,670 @@ +#include "DD4hep/DetFactoryHelper.h" +#include "TMatrixT.h" + +// todo: remove gaudi logging and properly capture output +#define endmsg std::endl +#define lLog std::cout +namespace MSG { +const std::string ERROR = " Error: "; +const std::string DEBUG = " Debug: "; +const std::string INFO = " Info: "; +} + +namespace det { + + unsigned ECalEndCapElementCounter = 0; + + unsigned ECalEndcapNumCalibLayers; + + double tForArcLength(double s, double bladeangle, double delZ, double r) { + + // some intermediate constants + double zpos = delZ/2.; + double zp = zpos/TMath::Tan(bladeangle); + double b = zp/(TMath::Sqrt(r*r-zp*zp)); + double c = (TMath::Tan(s/r) +b)/(1.-b*TMath::Tan(s/r)); + double d = c*c*r*r/(1+c*c); + return (TMath::Sqrt(d)-zp)*TMath::Sin(bladeangle); + + // try approximating the arclength as dx. Less accurate, but that + // approximation is used in calculating the LAr gap, so maybe this + // will make it more consistent? + //return s*TMath::Sin(bladeangle); + + } + + // return position of the inner edge of a blade + double getZmin(double r, double bladeangle, double delZ) { + // r: distance from the beamline + // bladeangle: angle of turbine blades wrt xy plane, in radians + // delZ: z extent of the blades + return TMath::Sqrt(r*r - ((delZ/2)/TMath::Tan(bladeangle))*((delZ/2)/TMath::Tan(bladeangle))); + } + + dd4hep::Solid buildOneBlade(double thickness_inner, + double thickness_outer, + double width, + double ro, double ri, + double bladeangle, + double delZ) + { + + dd4hep::Solid shapeBeforeSubtraction; + + // set max and min extent of the blade (along the z axis in the body frame) + double zmax = ro; + double zmin = getZmin(ri, bladeangle, delZ); + + dd4hep::Trd2 tmp1(thickness_inner/2., thickness_outer/2., width/2., width/2., (zmax-zmin)/2. ); + shapeBeforeSubtraction = tmp1; + + dd4hep::Tube allowedTube(ri, ro, delZ); + + return dd4hep::IntersectionSolid (shapeBeforeSubtraction, allowedTube, dd4hep::Transform3D(dd4hep::RotationZYX( 0, TMath::Pi()/2.-bladeangle, TMath::Pi()/2.),dd4hep::Position(0,0, -(zmin+zmax)/2.))); + + } + + void buildWheel(dd4hep::Detector& aLcdd, + dd4hep::SensitiveDetector& aSensDet, + dd4hep::Volume& aEnvelope, + dd4hep::xml::Handle_t& aXmlElement, + dd4hep::DetElement& bathDetElem, + float ri, float ro, float delZ, + unsigned iWheel) { + + + dd4hep::xml::DetElement calorimeterElem = aXmlElement.child(_Unicode(calorimeter)); + dd4hep::xml::DetElement genericBladeElem = calorimeterElem.child(_Unicode(turbineBlade)); + dd4hep::xml::DetElement absBladeElem = genericBladeElem.child(_Unicode(absorberBlade)); + dd4hep::xml::DetElement claddingElem = genericBladeElem.child(_Unicode(cladding)); + dd4hep::xml::DetElement glueElem = genericBladeElem.child(_Unicode(glue)); + dd4hep::xml::DetElement electrodeBladeElem = genericBladeElem.child(_Unicode(electrodeBlade)); + dd4hep::xml::DetElement nobleLiquidElem = genericBladeElem.child(_Unicode(nobleLiquidGap)); + + float BladeAngle = genericBladeElem.attr(_Unicode(angle)); + bool decreaseAnglePerWheel = genericBladeElem.attr(_Unicode(decreaseAnglePerWheel)); + lLog << MSG::DEBUG << "Making wheel with inner, outer radii " << ri << ", " << ro << std:: endl; + lLog << MSG::DEBUG << "Blade angle is " << BladeAngle << "; decrease angle per wheel? " << decreaseAnglePerWheel << endmsg; + dd4hep::xml::Dimension dim(aXmlElement.child(_Unicode(dimensions))); + double grmin = dim.rmin1(); + lLog << MSG::DEBUG << "delZ is " << delZ << endmsg; + if (decreaseAnglePerWheel) { + float tubeFracCovered = delZ/(2*grmin*TMath::Tan(BladeAngle)); + BladeAngle = TMath::ATan(delZ/(2*ri*tubeFracCovered)); + } + + if (TMath::Abs(TMath::Tan(BladeAngle)) < delZ/(2.*ri)) { + lLog << MSG::ERROR << "The requested blade angle is too small for the given delZ and ri values. Please adjust to at least " << TMath::ATan(delZ/(2.*ri))*180./TMath::Pi() << " degrees!" << endmsg; + return; + } + + Float_t xRange = delZ/(TMath::Sin(BladeAngle)); + + double delrPhiNoGap; + + float GlueThick = glueElem.attr(_Unicode(thickness)); + float CladdingThick = claddingElem.attr(_Unicode(thickness)); + float AbsThickMin = absBladeElem.attr(_Unicode(thickness))-(GlueThick+CladdingThick); + if (AbsThickMin < 0.) { + lLog << MSG::ERROR << "Error: requested absorber thickness is negative after accounting for glue and cladding thickness" << endmsg; + } + float ElectrodeThick = electrodeBladeElem.attr(_Unicode(thickness)); + float LArgapi = nobleLiquidElem.attr(_Unicode(gap)); + + bool sameNUnitCells = genericBladeElem.attr(_Unicode(sameNUnitCells)); + auto nUnitCellsStrArr = genericBladeElem.attr(_Unicode(nUnitCells)); + char* nUnitCellsCStr = strtok(const_cast(nUnitCellsStrArr.c_str()), " "); + int nUnitCells = 0; + if (!sameNUnitCells) { + for (unsigned i = 0; i < iWheel; i++) { + nUnitCellsCStr = strtok(NULL, " "); + } + std::string nUnitCellsStr = nUnitCellsCStr; + nUnitCells = std::stoi(nUnitCellsStr); + } + int nUnitCellsLeastCommonMultiple = genericBladeElem.attr(_Unicode(nUnitCellsLeastCommonMultiple)); + + bool scaleBladeThickness = absBladeElem.attr(_Unicode(scaleThickness)); + float bladeThicknessScaleFactor = absBladeElem.attr(_Unicode(thicknessScaleFactor)); + + lLog << MSG::DEBUG << "nUnitCells: " << nUnitCells << endmsg; + + float AbsThicki = AbsThickMin; + // make volumes for the noble liquid, electrode, and absorber blades + float AbsThicko; + if (scaleBladeThickness) { + AbsThicko = AbsThicki + bladeThicknessScaleFactor*((ro/ri)-1.)*AbsThicki; + } else { + AbsThicko = AbsThicki; + } + + // Calculate gap thickness at inner layer + double circ = 2*TMath::Pi()*ri; + double x2 =(AbsThickMin+(GlueThick+CladdingThick)+ElectrodeThick)/TMath::Sin(BladeAngle); + double y2 = TMath::Sqrt(ri*ri-x2*x2); + double rPhi1 = ri*TMath::Pi()/2.; + double rPhi2 = ri*TMath::ATan(y2/x2); + delrPhiNoGap = TMath::Abs(rPhi1-rPhi2); + double leftoverS = (circ - nUnitCells*delrPhiNoGap); + double delrPhiGapOnly = leftoverS/(2*nUnitCells); + LArgapi = delrPhiGapOnly*TMath::Sin(BladeAngle); + lLog << MSG::DEBUG << "LArGap at inner radius is " << LArgapi << endmsg; + + // now find gap at outer radius + circ = 2*TMath::Pi()*ro; + x2 = (AbsThicko+GlueThick+CladdingThick+ElectrodeThick)/TMath::Sin(BladeAngle); + y2 = TMath::Sqrt(ro*ro-x2*x2); + rPhi1 = ro*TMath::Pi()/2.; + rPhi2 = ro*TMath::ATan(y2/x2); + delrPhiNoGap = TMath::Abs(rPhi1-rPhi2); + leftoverS = (circ - nUnitCells*delrPhiNoGap); + delrPhiGapOnly = leftoverS/(2*nUnitCells); + float LArgapo = delrPhiGapOnly*TMath::Sin(BladeAngle); + // LArgapo *= 2.; + + float riLayer = ri; + + std::vector claddingLayerVols; + std::vector glueLayerVols; + std::vector absBladeLayerVols; + std::vector LArTotalLayerVols; + std::vector electrodeBladeLayerVols; + + + dd4hep::Solid passiveShape = buildOneBlade(AbsThicki+GlueThick+CladdingThick, AbsThicko+GlueThick+CladdingThick, xRange, ro, ri, BladeAngle, delZ ); + dd4hep::Volume passiveVol("passive", passiveShape, aLcdd.material("Air")); + + dd4hep::Solid activeShape = buildOneBlade(ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2, xRange, ro, ri, BladeAngle, delZ); + dd4hep::Volume activeVol("active", activeShape, aLcdd.material("Air")); + + unsigned numNonActiveLayers = 1; + // check that either all non-active volumes are set to sensitive (for + // sampling fraction calculations) or none are (for normal running) + bool allNonActiveSensitive = ( claddingElem.isSensitive() && + glueElem.isSensitive() && + absBladeElem.isSensitive() && + electrodeBladeElem.isSensitive() ); + bool allNonActiveNotSensitive = ( !claddingElem.isSensitive() && + !glueElem.isSensitive() && + !absBladeElem.isSensitive() && + !electrodeBladeElem.isSensitive() ); + if (allNonActiveSensitive) { + numNonActiveLayers = ECalEndcapNumCalibLayers; + } else if (allNonActiveNotSensitive) { + numNonActiveLayers = 1; + } else { + lLog << MSG::ERROR << "Some non-active layers are sensitive and others are not -- this is likely a misconfiguration"; + } + + float delrNonActive = (ro-ri)/numNonActiveLayers; + float delrActive = (ro-ri)/ECalEndcapNumCalibLayers; + + for (unsigned iLayer = 0; iLayer < numNonActiveLayers; iLayer++) { + float roLayer = riLayer + delrNonActive; + lLog << MSG::INFO << "Making layer in inner, outer radii " << riLayer << " " << roLayer << endmsg; + + if (scaleBladeThickness) { + AbsThicko = AbsThicki + bladeThicknessScaleFactor*((roLayer/riLayer)-1.)*AbsThicki; + } else { + AbsThicko = AbsThicki; + } + lLog << MSG::DEBUG << "Inner and outer absorber thicknesses " << AbsThicki << " " << AbsThicko << endmsg; + dd4hep::Solid claddingLayer = buildOneBlade(AbsThicki+GlueThick+CladdingThick, AbsThicko+GlueThick+CladdingThick, xRange, roLayer, riLayer, BladeAngle, delZ ); + + dd4hep::Solid glueLayer = buildOneBlade(AbsThicki+GlueThick, AbsThicko+GlueThick, xRange, roLayer, riLayer, BladeAngle, delZ ); + + // dd4hep::SubtractionSolid claddingLayer(absGlueCladdingLayer, absGlueLayer); + dd4hep::Solid absBladeLayer = buildOneBlade(AbsThicki, AbsThicko, xRange, roLayer, riLayer, BladeAngle, delZ ); + + // dd4hep::SubtractionSolid glueLayer(absGlueLayer, absBladeLayer); + dd4hep::Volume claddingLayerVol("claddingLayer", claddingLayer, aLcdd.material(claddingElem.materialStr())); + if (claddingElem.isSensitive()) { + claddingLayerVol.setSensitiveDetector(aSensDet); + } + claddingLayerVols.push_back(claddingLayerVol); + + dd4hep::Volume glueLayerVol("glueLayer", glueLayer, aLcdd.material(glueElem.materialStr())); + if (glueElem.isSensitive()) { + glueLayerVol.setSensitiveDetector(aSensDet); + } + glueLayerVols.push_back(glueLayerVol); + + dd4hep::Volume absBladeLayerVol("absBladeLayer", absBladeLayer, aLcdd.material(absBladeElem.materialStr())); + if (absBladeElem.isSensitive()) { + absBladeLayerVol.setSensitiveDetector(aSensDet); + } + absBladeLayerVols.push_back(absBladeLayerVol); + + riLayer = roLayer; + AbsThicki = AbsThicko; + } + + riLayer = ri; + + AbsThicki = AbsThickMin; + + for (unsigned iLayer = 0; iLayer < ECalEndcapNumCalibLayers; iLayer++) { + + float roLayer = riLayer + delrActive; + + if (scaleBladeThickness) { + AbsThicko = AbsThicki + bladeThicknessScaleFactor*((roLayer/riLayer)-1.)*AbsThicki; + } else { + AbsThicko = AbsThicki; + } + + + // now find gap at outer layer + circ = 2*TMath::Pi()*roLayer; + x2 = (AbsThicko+GlueThick+CladdingThick+ElectrodeThick)/TMath::Sin(BladeAngle); + y2 = TMath::Sqrt(roLayer*roLayer-x2*x2); + rPhi1 = roLayer*TMath::Pi()/2.; + rPhi2 = roLayer*TMath::ATan(y2/x2); + delrPhiNoGap = TMath::Abs(rPhi1-rPhi2); + leftoverS = (circ - nUnitCells*delrPhiNoGap); + delrPhiGapOnly = leftoverS/(2*nUnitCells); + LArgapo = delrPhiGapOnly*TMath::Sin(BladeAngle); + lLog << MSG::DEBUG << "Outer LAr gap is " << LArgapo << endmsg ; + lLog << MSG::INFO << "Inner and outer thicknesses of noble liquid volume " << ElectrodeThick+LArgapi*2 << " " << ElectrodeThick+LArgapo*2 << endmsg; + + dd4hep::Solid electrodeBladeAndGapLayer = buildOneBlade(ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2, xRange, roLayer, riLayer, BladeAngle, delZ); + + dd4hep::Solid electrodeBladeLayer = buildOneBlade(ElectrodeThick, ElectrodeThick, xRange, roLayer, riLayer, BladeAngle, delZ); + + dd4hep::Volume electrodeBladeLayerVol("electrodeBladeLayer", electrodeBladeLayer, aLcdd.material(electrodeBladeElem.materialStr())); + if (electrodeBladeElem.isSensitive()) { + electrodeBladeLayerVol.setSensitiveDetector(aSensDet); + } + electrodeBladeLayerVols.push_back(electrodeBladeLayerVol); + + // dd4hep::SubtractionSolid LArShapeTotalLayer(electrodeBladeAndGapLayer, electrodeBladeLayer); + dd4hep::Volume LArTotalLayerVol("LArTotalLayerVol", electrodeBladeAndGapLayer, aLcdd.material(nobleLiquidElem.materialStr())); + + if ( nobleLiquidElem.isSensitive() ) { + LArTotalLayerVol.setSensitiveDetector(aSensDet); + } + LArTotalLayerVols.push_back(LArTotalLayerVol); + + riLayer = roLayer; + LArgapi = LArgapo; + AbsThicki = AbsThicko; + } + lLog << MSG::INFO << "ECal endcap materials: nobleLiquid: " << nobleLiquidElem.materialStr() << " absorber: " << absBladeElem.materialStr() << " electrode: " << electrodeBladeElem.materialStr() << endmsg; + + int nUnitCellsToDraw = nUnitCells; + // nUnitCellsToDraw = 2; + + lLog << MSG::INFO << "Number of unit cells "<< nUnitCells << endmsg; + + // place all components of the absorber blade inside passive volume + + unsigned iLayer = 0; + + riLayer = ri; + + for (auto absBladeLayerVol: absBladeLayerVols) { + + float roLayer = riLayer+delrNonActive; + + dd4hep::Position posLayer(0,0,(riLayer-ri+roLayer-ro)/2.); + dd4hep::PlacedVolume absBladeVol_pv = glueLayerVols[iLayer].placeVolume(absBladeLayerVol, posLayer); + + absBladeVol_pv.addPhysVolID("subtype", 0); // 0 = absorber, 1 = glue, 2 = cladding + lLog << MSG::DEBUG << "Blade layer, rho is " << iLayer << " " << absBladeVol_pv.position().Rho() << " " << roLayer/2. << endmsg; + absBladeVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + riLayer = roLayer; + iLayer++; + } + + riLayer = ri; + iLayer =0; + + for (auto glueLayerVol: glueLayerVols) { + + float roLayer = riLayer+delrNonActive; + + dd4hep::Position posLayer(0,0,(riLayer-ri+roLayer-ro)/2.); + dd4hep::PlacedVolume glueVol_pv = claddingLayerVols[iLayer].placeVolume(glueLayerVol, posLayer); + + + glueVol_pv.addPhysVolID("subtype", 1); // 0 = absorber, 1 = glue, 2 = cladding + glueVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + // dd4hep::DetElement glueDetElem(passiveDetElem, "glue_",ECalEndCapElementCounter++); + // glueDetElem.setPlacement(glueVol_pv); + + riLayer = roLayer; + iLayer++; + } + + riLayer = ri; + iLayer =0; + + double zminri = getZmin(ri, BladeAngle, delZ); + + for (auto claddingLayerVol: claddingLayerVols) { + + float roLayer = riLayer+delrNonActive; + + double zminLayer = getZmin(riLayer, BladeAngle, delZ); + + dd4hep::Position posLayer(0,0,(zminLayer-zminri+roLayer-ro)/2.); + dd4hep::PlacedVolume claddingVol_pv = passiveVol.placeVolume(claddingLayerVol, posLayer); + + claddingVol_pv.addPhysVolID("subtype", 2); // 0 = absorber, 1 = glue, 2 = cladding + claddingVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + // dd4hep::DetElement claddingDetElem(passiveDetElem, "cladding_", ECalEndCapElementCounter++); + // claddingDetElem.setPlacement(claddingVol_pv); + + riLayer = roLayer; + iLayer++; + } + + + riLayer = ri; + iLayer = 0; + + for (auto electrodeBladeLayerVol: electrodeBladeLayerVols) { + + float roLayer = riLayer+delrActive; + + dd4hep::PlacedVolume electrodeBladeVol_pv = LArTotalLayerVols[iLayer].placeVolume(electrodeBladeLayerVol); + electrodeBladeVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + riLayer = roLayer; + iLayer++; + } + + riLayer = ri; + iLayer = 0; + + for (auto LArTotalLayerVol: LArTotalLayerVols) { + + float roLayer = riLayer+delrActive; + + double zminLayer = getZmin(riLayer, BladeAngle, delZ); + + dd4hep::Position posLayer(0,0,(zminLayer-zminri+roLayer-ro)/2.); + std::cout << "for active, riLayer, ri, roLayer, ro = " << riLayer << " " << ri << " " << roLayer << " " << ro << std::endl; + + dd4hep::PlacedVolume LArVol_pv(activeVol.placeVolume(LArTotalLayerVol, posLayer)); + lLog << MSG::DEBUG << "LAr layer: " << iLayer << endmsg; + LArVol_pv.addPhysVolID("layer", iWheel*ECalEndcapNumCalibLayers+iLayer); + + riLayer = roLayer; + iLayer++; + } + + for (int iUnitCell = 0; iUnitCell < nUnitCellsToDraw; iUnitCell++) { + + int modIndex = iUnitCell-nUnitCellsToDraw/2; + if (modIndex < 0) modIndex += nUnitCells; + float phi = (iUnitCell-nUnitCellsToDraw/2)*2*TMath::Pi()/nUnitCells; + float delPhi = 2*TMath::Pi()/nUnitCells; + + lLog << MSG::DEBUG << "Placing blade, ro, ri = " << ro << " " << ri << endmsg; + TGeoRotation tgr; + tgr.RotateZ(BladeAngle*180/TMath::Pi()); + tgr.RotateX(-phi*180/TMath::Pi()); + tgr.RotateY(90); + + const Double_t *rotMatPtr; + + rotMatPtr = tgr.GetRotationMatrix(); + TMatrixT rotMat(3,3, rotMatPtr); + dd4hep::Rotation3D r3d; + r3d.SetComponents(rotMat(0,0), rotMat(0,1), rotMat(0,2), + rotMat(1,0), rotMat(1,1), rotMat(1,2), + rotMat(2,0), rotMat(2,1), rotMat(2,2)); + + tgr.Clear(); + tgr.RotateZ(BladeAngle*180/TMath::Pi()); + tgr.RotateX(-(phi+delPhi/2.)*180/TMath::Pi()); + tgr.RotateY(90); + + rotMatPtr = tgr.GetRotationMatrix(); + TMatrixT rotMat2(3,3, rotMatPtr); + dd4hep::Rotation3D r3d2; + r3d2.SetComponents(rotMat2(0,0), rotMat2(0,1), rotMat2(0,2), + rotMat2(1,0), rotMat2(1,1), rotMat2(1,2), + rotMat2(2,0), rotMat2(2,1), rotMat2(2,2)); + + riLayer = ri; + + float xCell = ((ro+zminri)/2.)*TMath::Cos(phi); + float yCell = ((ro+zminri)/2.)*TMath::Sin(phi); //ri*TMath::Sin(phi)/6.; + float zCell = 0.; + + dd4hep::Transform3D comCell(r3d, dd4hep::Translation3D(xCell,yCell,zCell)); + + // place passive volume in LAr bath + dd4hep::PlacedVolume passivePhysVol = aEnvelope.placeVolume(passiveVol, comCell); + passivePhysVol.addPhysVolID("module", modIndex*nUnitCellsLeastCommonMultiple/nUnitCells); + passivePhysVol.addPhysVolID("wheel", iWheel); + passivePhysVol.addPhysVolID("type", 1); // 0 = active, 1 = passive, 2 = readout + dd4hep::DetElement passiveDetElem(bathDetElem, "passive_" + std::to_string(iUnitCell)+"_"+std::to_string(iWheel), ECalEndCapElementCounter++); + passiveDetElem.setPlacement(passivePhysVol); + + // place active volume in LAr bath + + xCell = ((ro+zminri)/2.)*TMath::Cos(phi+delPhi/2.); + yCell = ((ro+zminri)/2.)*TMath::Sin(phi+delPhi/2.); //ri*TMath::Sin(phi)/6.; + zCell = 0.; + dd4hep::Transform3D comCell2(r3d2, dd4hep::Translation3D(xCell,yCell,zCell)); + dd4hep::PlacedVolume activePhysVol = aEnvelope.placeVolume(activeVol, comCell2); + activePhysVol.addPhysVolID("module", modIndex*nUnitCellsLeastCommonMultiple/nUnitCells); + activePhysVol.addPhysVolID("wheel", iWheel); + activePhysVol.addPhysVolID("type", 0); // 0 = active, 1 = passive, 2 = readout + + dd4hep::DetElement activeDetElem(bathDetElem, "active_" + std::to_string(iUnitCell)+"_"+std::to_string(iWheel), ECalEndCapElementCounter++); + activeDetElem.setPlacement(activePhysVol); + + riLayer = ri; + iLayer =0; + + + + lLog << MSG::DEBUG << "LArTotalLayerVols.size = " << LArTotalLayerVols.size() << endmsg; + + } + + return; + + } + + void buildOneSide_Turbine(dd4hep::Detector& aLcdd, dd4hep::SensitiveDetector& aSensDet, + dd4hep::Volume& aEnvelope, dd4hep::xml::Handle_t& aXmlElement, + unsigned& iModule) { + + dd4hep::xml::DetElement calo = aXmlElement.child(_Unicode(calorimeter)); + dd4hep::xml::Dimension caloDim(calo.dimensions()); + + + dd4hep::xml::DetElement blade = calo.child(_Unicode(turbineBlade)); + dd4hep::xml::DetElement nobleLiquid = blade.child(_Unicode(nobleLiquidGap)); + + dd4hep::xml::DetElement xmlDetElem = aXmlElement; + std::string nameDet = xmlDetElem.nameStr(); + dd4hep::DetElement caloDetElem(nameDet, xmlDetElem.id()); + + dd4hep::xml::Dimension dim(aXmlElement.child(_Unicode(dimensions))); + + //build cryostat +// Retrieve cryostat data + dd4hep::xml::DetElement cryostat = calo.child(_Unicode(cryostat)); + dd4hep::xml::Dimension cryoDim(cryostat.dimensions()); + double cryoThicknessFront = cryoDim.rmin2() - cryoDim.rmin1(); + + dd4hep::xml::DetElement cryoFront = cryostat.child(_Unicode(front)); + dd4hep::xml::DetElement cryoBack = cryostat.child(_Unicode(back)); + dd4hep::xml::DetElement cryoSide = cryostat.child(_Unicode(side)); + bool cryoFrontSensitive = cryoFront.isSensitive(); + bool cryoBackSensitive = cryoBack.isSensitive(); + bool cryoSideSensitive = cryoSide.isSensitive(); + + double bathRmin = caloDim.rmin(); // - margin for inclination + double bathRmax = caloDim.rmax(); // + margin for inclination + double bathDelZ = caloDim.dz(); + dd4hep::Tube bathOuterShape(bathRmin, bathRmax, bathDelZ); // make it 4 volumes + 5th for detector envelope + dd4hep::Tube bathAndServicesOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), caloDim.dz()); // make it 4 volumes + 5th for detector envelope + + lLog << MSG::INFO << "Cryostat front thickness is " << cryoDim.rmin2() << endmsg; + if (cryoThicknessFront > 0) { + // 1. Create cryostat + dd4hep::Tube cryoFrontShape(cryoDim.rmin1(), cryoDim.rmin2(), cryoDim.dz()); + dd4hep::Tube cryoBackShape(cryoDim.rmax1(), cryoDim.rmax2(), cryoDim.dz()); + dd4hep::Tube cryoSideOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz()); + dd4hep::SubtractionSolid cryoSideShape(cryoSideOuterShape, bathAndServicesOuterShape); + lLog << MSG::INFO << "ECAL endcap cryostat: front: rmin (cm) = " << cryoDim.rmin1() << " rmax (cm) = " << cryoDim.rmin2() << " dz (cm) = " << cryoDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL encdap cryostat: back: rmin (cm) = " << cryoDim.rmax1() << " rmax (cm) = " << cryoDim.rmax2() << " dz (cm) = " << cryoDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL endcap cryostat: side: rmin (cm) = " << cryoDim.rmin2() << " rmax (cm) = " << cryoDim.rmax1() << " dz (cm) = " << cryoDim.dz() - caloDim.dz() << endmsg; + lLog << MSG::INFO << "Cryostat is made out of " << cryostat.materialStr() << endmsg; + + dd4hep::Volume cryoFrontVol(cryostat.nameStr()+"_front", cryoFrontShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoBackVol(cryostat.nameStr()+"_back", cryoBackShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoSideVol(cryostat.nameStr()+"_side", cryoSideShape, aLcdd.material(cryostat.materialStr())); + dd4hep::PlacedVolume cryoFrontPhysVol = aEnvelope.placeVolume(cryoFrontVol); + dd4hep::PlacedVolume cryoBackPhysVol = aEnvelope.placeVolume(cryoBackVol); + dd4hep::PlacedVolume cryoSidePhysVol = aEnvelope.placeVolume(cryoSideVol); + unsigned sidetype = 0x4; // probably not needed anymore... + if (cryoFrontSensitive) { + cryoFrontVol.setSensitiveDetector(aSensDet); + cryoFrontPhysVol.addPhysVolID("cryo", 1); + cryoFrontPhysVol.addPhysVolID("type", sidetype+1); + lLog << MSG::INFO << "Cryostat front volume set as sensitive" << endmsg; + } + if (cryoBackSensitive) { + cryoBackVol.setSensitiveDetector(aSensDet); + cryoBackPhysVol.addPhysVolID("cryo", 1); + cryoBackPhysVol.addPhysVolID("type", sidetype+2); + lLog << MSG::INFO << "Cryostat back volume set as sensitive" << endmsg; + } + if (cryoSideSensitive) { + cryoSideVol.setSensitiveDetector(aSensDet); + cryoSidePhysVol.addPhysVolID("cryo", 1); + cryoSidePhysVol.addPhysVolID("type", sidetype+3); + lLog << MSG::INFO << "Cryostat front volume set as sensitive" << endmsg; + } + dd4hep::DetElement cryoFrontDetElem(caloDetElem, "cryo_front", 0); + cryoFrontDetElem.setPlacement(cryoFrontPhysVol); + dd4hep::DetElement cryoBackDetElem(caloDetElem, "cryo_back", 0); + cryoBackDetElem.setPlacement(cryoBackPhysVol); + dd4hep::DetElement cryoSideDetElem(caloDetElem, "cryo_side", 0); + cryoSideDetElem.setPlacement(cryoSidePhysVol); + } + + // 2. Create noble liquid bath + std::string nobleLiquidMaterial = nobleLiquid.materialStr(); + dd4hep::Volume bathVol(nobleLiquidMaterial + "_bath", bathOuterShape, aLcdd.material(nobleLiquidMaterial)); + lLog << MSG::INFO << "ECAL endcap bath: material = " << nobleLiquidMaterial << " rmin (cm) = " << bathRmin + << " rmax (cm) = " << bathRmax << " dz (cm) = " << caloDim.dz() << " thickness in front of ECal (cm) = " << caloDim.rmin() - cryoDim.rmin2() + << " thickness behind ECal (cm) = " << cryoDim.rmax1() - caloDim.rmax() << endmsg; + dd4hep::DetElement bathDetElem(caloDetElem, "bath", 1); + + // 3. Create detector structure + double length = dim.dz() * 2.; + double zOffsetEnvelope = -length / 2.; + + dd4hep::xml::DetElement supportTubeElem = calo.child(_Unicode(supportTube)); + unsigned nWheels = supportTubeElem.attr(_Unicode(nWheels)); + lLog << MSG::INFO << "Will build " << nWheels << " wheels" << endmsg; + double rmin = bathRmin; + double rmax = bathRmax; + float radiusRatio = pow(rmax/rmin, 1./nWheels); + double ro = rmin*radiusRatio; + double ri = rmin; + + float supportTubeThickness=supportTubeElem.thickness(); + + for (unsigned iWheel = 0; iWheel < nWheels; iWheel++) { + + dd4hep::Tube supportTube(ro, ro+supportTubeThickness, bathDelZ ); + + dd4hep::Volume supportTubeVol("supportTube", supportTube, aLcdd.material(supportTubeElem.materialStr())); + if (supportTubeElem.isSensitive()) { + supportTubeVol.setSensitiveDetector(aSensDet); + } + dd4hep::PlacedVolume supportTube_pv = bathVol.placeVolume(supportTubeVol, dd4hep::Position(0,0,zOffsetEnvelope + dim.dz() )); + supportTube_pv.addPhysVolID("cryo", 1); + // supportTube_pv.addPhysVolID("side",sign); + supportTube_pv.addPhysVolID("wheel", iWheel); + dd4hep::DetElement supportTubeDetElem(bathDetElem, "supportTube_"+std::to_string(iWheel), 0); + supportTubeDetElem.setPlacement(supportTube_pv); + + + buildWheel(aLcdd, aSensDet, bathVol, aXmlElement, bathDetElem, ri+supportTubeThickness, ro, bathDelZ*2, iWheel); + ri = ro; + ro *= radiusRatio; + if (ro > rmax) ro = rmax; + } + + + dd4hep::PlacedVolume bathPhysVol = aEnvelope.placeVolume(bathVol); + bathDetElem.setPlacement(bathPhysVol); + + lLog << MSG::DEBUG << "Total number of modules: " << iModule << endmsg; + + return; +} + + + + static dd4hep::Ref_t +createECalEndcapTurbine(dd4hep::Detector& aLcdd, dd4hep::xml::Handle_t aXmlElement, dd4hep::SensitiveDetector aSensDet) { + + dd4hep::xml::DetElement xmlDetElem = aXmlElement; + std::string nameDet = xmlDetElem.nameStr(); + int idDet = xmlDetElem.id(); + dd4hep::xml::Dimension dim(xmlDetElem.dimensions()); + dd4hep::DetElement caloDetElem(nameDet, idDet); + dd4hep::xml::Dimension sdType = xmlDetElem.child(_U(sensitive)); + aSensDet.setType(sdType.typeStr()); + + ECalEndcapNumCalibLayers = aLcdd.constant("ECalEndcapNumCalibLayers"); + + + // Create air envelope for one endcap (will be copied to make both endcaps) + dd4hep::Tube endcapShape( dim.rmin1(), dim.rmax1(), dim.dz()); + + dd4hep::Volume envelopeVol(nameDet + "_vol", endcapShape, aLcdd.material("Air")); + + + // dd4hep::DetElement caloPositiveDetElem(caloDetElem, "positive", 0); + // dd4hep::DetElement caloNegativeDetElem(caloDetElem, "negative", 0); + + lLog << MSG::DEBUG << "Placing dector on the positive side: (cm) " << dim.z_offset() << " with min, max radii " << dim.rmin1() << " " << dim.rmax1() << endmsg; + unsigned iModule = 0; + buildOneSide_Turbine(aLcdd, aSensDet, envelopeVol, aXmlElement, iModule); + // lLog << MSG::DEBUG << "Placing dector on the negative side: (cm) " << -dim.z_offset() << " with min, max radii " << dim.rmin1() << " " << dim.rmax() << endmsg; + // buildOneSide_Turbine(aLcdd, aSensDet, envelopeNegativeVol, aXmlElement, -1, iModule); + + dd4hep::Assembly endcapsAssembly("ECalEndcaps_turbine"); + + // Place the envelope + dd4hep::Transform3D envelopePositiveVolume_tr(dd4hep::RotationZYX( 0 ,0,0), dd4hep::Translation3D(0, 0, dim.z_offset())); + dd4hep::PlacedVolume envelopePositivePhysVol = endcapsAssembly.placeVolume(envelopeVol, envelopePositiveVolume_tr); + envelopePositivePhysVol.addPhysVolID("side", 1); + + dd4hep::DetElement caloPositiveDetElem(caloDetElem, "positive", 0); + caloPositiveDetElem.setPlacement(envelopePositivePhysVol); + + // make another placement for the negative z endcap + dd4hep::Transform3D envelopeNegativeVolume_tr(dd4hep::RotationZYX( 0 ,0,180*dd4hep::deg), dd4hep::Translation3D(0, 0, -dim.z_offset())); + dd4hep::PlacedVolume envelopeNegativePhysVol = + endcapsAssembly.placeVolume(envelopeVol, envelopeNegativeVolume_tr); + envelopeNegativePhysVol.addPhysVolID("side", -1); + + dd4hep::DetElement caloNegativeDetElem(caloDetElem, "negative", 0); + caloNegativeDetElem.setPlacement(envelopeNegativePhysVol); + + dd4hep::Volume motherVol = aLcdd.pickMotherVolume(caloDetElem); + dd4hep::PlacedVolume envelopePhysVol = motherVol.placeVolume(endcapsAssembly); + caloDetElem.setPlacement(envelopePhysVol); + envelopePhysVol.addPhysVolID("system", idDet); + return caloDetElem; +} +} // namespace det + +DECLARE_DETELEMENT(ECalEndcap_Turbine_o1_v01, det::createECalEndcapTurbine) diff --git a/detector/calorimeter/ECalEndcap_o1_v01_geo.cpp b/detector/calorimeter/ECalEndcap_o1_v01_geo.cpp index f27941378..b08119890 100644 --- a/detector/calorimeter/ECalEndcap_o1_v01_geo.cpp +++ b/detector/calorimeter/ECalEndcap_o1_v01_geo.cpp @@ -12,9 +12,7 @@ //========================================================================= #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "TMath.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/ECalPlug_o1_v01_geo.cpp b/detector/calorimeter/ECalPlug_o1_v01_geo.cpp index be6b57242..ea9791421 100644 --- a/detector/calorimeter/ECalPlug_o1_v01_geo.cpp +++ b/detector/calorimeter/ECalPlug_o1_v01_geo.cpp @@ -12,9 +12,6 @@ //========================================================================= #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" -#include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" using namespace std; diff --git a/detector/calorimeter/Hcal_BarrelSD_v00.cpp b/detector/calorimeter/Hcal_BarrelSD_v00.cpp index 8ee66e3e2..d64b308cc 100644 --- a/detector/calorimeter/Hcal_BarrelSD_v00.cpp +++ b/detector/calorimeter/Hcal_BarrelSD_v00.cpp @@ -8,10 +8,7 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "XML/Layering.h" #include "XML/Utilities.h" -#include "DDRec/DetectorData.h" -#include "LcgeoExceptions.h" using namespace std; diff --git a/detector/calorimeter/Hcal_Barrel_SD_v01.cpp b/detector/calorimeter/Hcal_Barrel_SD_v01.cpp index 8da93799f..feb38eb70 100644 --- a/detector/calorimeter/Hcal_Barrel_SD_v01.cpp +++ b/detector/calorimeter/Hcal_Barrel_SD_v01.cpp @@ -5,10 +5,8 @@ // $Id: $ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "XML/Layering.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" -#include "LcgeoExceptions.h" using namespace std; diff --git a/detector/calorimeter/Hcal_Barrel_SD_v02.cpp b/detector/calorimeter/Hcal_Barrel_SD_v02.cpp index f186dba6c..acc3ff518 100644 --- a/detector/calorimeter/Hcal_Barrel_SD_v02.cpp +++ b/detector/calorimeter/Hcal_Barrel_SD_v02.cpp @@ -5,10 +5,8 @@ // $Id: $ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "XML/Layering.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" -#include "LcgeoExceptions.h" using namespace std; diff --git a/detector/calorimeter/Hcal_EndcapRing_SD_v01.cpp b/detector/calorimeter/Hcal_EndcapRing_SD_v01.cpp index 97ade488f..620ad7ac6 100644 --- a/detector/calorimeter/Hcal_EndcapRing_SD_v01.cpp +++ b/detector/calorimeter/Hcal_EndcapRing_SD_v01.cpp @@ -41,7 +41,6 @@ */ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "DD4hep/Shapes.h" #include "XML/Utilities.h" diff --git a/detector/calorimeter/Hcal_Endcaps_SD_v01.cpp b/detector/calorimeter/Hcal_Endcaps_SD_v01.cpp index 9243c7555..203298ae7 100644 --- a/detector/calorimeter/Hcal_Endcaps_SD_v01.cpp +++ b/detector/calorimeter/Hcal_Endcaps_SD_v01.cpp @@ -42,7 +42,6 @@ */ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "DD4hep/Shapes.h" #include "XML/Utilities.h" diff --git a/detector/calorimeter/Hcal_Endcaps_SD_v02.cpp b/detector/calorimeter/Hcal_Endcaps_SD_v02.cpp index 36bf73289..b9633d400 100644 --- a/detector/calorimeter/Hcal_Endcaps_SD_v02.cpp +++ b/detector/calorimeter/Hcal_Endcaps_SD_v02.cpp @@ -34,7 +34,6 @@ // ² //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/README.md b/detector/calorimeter/README.md index 3bcd9a375..ddb6a0b1b 100644 --- a/detector/calorimeter/README.md +++ b/detector/calorimeter/README.md @@ -22,6 +22,13 @@ This sub-detector makes calorimeter endcaps (original and reflected). It is used ### o1_v01 Original version taken from [FCCDetectors](https://github.com/HEP-FCC/FCCDetectors/blob/70a989a6fc333610e3b1b979c3596da9c41543d8/Detector/DetFCChhCalDiscs/src/CaloEndcapDiscs_geo.cpp). +## ECalEndcap_Turbine + +Sub-detector for the ecal endcaps, with the absorbers and readout boards arranged in a "turbine-like" geometry. + +### 01_v01 +Initial implementation. A custom segmentation that creates readout cells and constant radius and z is also included (FCCSWEndcapTurbine_k4geo). + ## HCalTileBarrel This sub-detector makes calorimeter barrel. It is used in ALLEGRO detector concept. @@ -35,4 +42,8 @@ This sub-detector makes calorimeter endcaps. Each endcap is made up by three cyl Original version taken from [FCCDetectors](https://github.com/HEP-FCC/FCCDetectors/blob/70a989a6fc333610e3b1b979c3596da9c41543d8/Detector/DetFCCeeHCalTile/src/HCalThreePartsEndcap_geo.cpp#L4). +## dual-readout +### o1_v01 +This sub-detector makes full 4-pi monolithic fiber dual-readout calorimeter. +Inside the single tower (trapezoidal copper absorber), two types of optical fibers (Cherenkov and scintillation) are implemented. The readout (SiPM) is attached at the rear side of the tower. The tower is repeated in both eta and phi direction to cover both barrel and endcap region. diff --git a/detector/calorimeter/SECalEndcap_o1_v01_geo.cpp b/detector/calorimeter/SECalEndcap_o1_v01_geo.cpp index bf16d0348..8b89582c9 100644 --- a/detector/calorimeter/SECalEndcap_o1_v01_geo.cpp +++ b/detector/calorimeter/SECalEndcap_o1_v01_geo.cpp @@ -10,10 +10,8 @@ //========================================================================= #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "DD4hep/DetType.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "TMath.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/SEcal04_Barrel.cpp b/detector/calorimeter/SEcal04_Barrel.cpp index 0470d9fe7..7f139de1f 100644 --- a/detector/calorimeter/SEcal04_Barrel.cpp +++ b/detector/calorimeter/SEcal04_Barrel.cpp @@ -9,9 +9,7 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/SEcal04_Barrel_v01.cpp b/detector/calorimeter/SEcal04_Barrel_v01.cpp index e97aa5565..5f5ad923b 100644 --- a/detector/calorimeter/SEcal04_Barrel_v01.cpp +++ b/detector/calorimeter/SEcal04_Barrel_v01.cpp @@ -9,9 +9,7 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/SEcal04_ECRing.cpp b/detector/calorimeter/SEcal04_ECRing.cpp index c5e8aa414..599ebf117 100644 --- a/detector/calorimeter/SEcal04_ECRing.cpp +++ b/detector/calorimeter/SEcal04_ECRing.cpp @@ -37,7 +37,6 @@ */ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/SEcal05_Barrel.cpp b/detector/calorimeter/SEcal05_Barrel.cpp index a0329b4e2..220e73f81 100644 --- a/detector/calorimeter/SEcal05_Barrel.cpp +++ b/detector/calorimeter/SEcal05_Barrel.cpp @@ -9,16 +9,12 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" -#include "DDSegmentation/MegatileLayerGridXY.h" -#include "DDSegmentation/WaferGridXY.h" #include "SEcal05_Helpers.h" diff --git a/detector/calorimeter/SEcal05_ECRing.cpp b/detector/calorimeter/SEcal05_ECRing.cpp index 186ad7605..69f53be00 100644 --- a/detector/calorimeter/SEcal05_ECRing.cpp +++ b/detector/calorimeter/SEcal05_ECRing.cpp @@ -73,7 +73,6 @@ Start ECRing ... #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/SEcal05_Endcaps.cpp b/detector/calorimeter/SEcal05_Endcaps.cpp index a0f4c69ab..ee4bce3e2 100644 --- a/detector/calorimeter/SEcal05_Endcaps.cpp +++ b/detector/calorimeter/SEcal05_Endcaps.cpp @@ -40,7 +40,6 @@ #include "DD4hep/Shapes.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" -#include "DDSegmentation/WaferGridXY.h" #include using namespace std; diff --git a/detector/calorimeter/SEcal05_Helpers.h b/detector/calorimeter/SEcal05_Helpers.h index 80e8e829a..bfea1f4fe 100644 --- a/detector/calorimeter/SEcal05_Helpers.h +++ b/detector/calorimeter/SEcal05_Helpers.h @@ -6,9 +6,7 @@ #include "DDRec/DetectorData.h" #include "XML/Layering.h" -#include "XML/Utilities.h" -#include "DD4hep/Segmentations.h" #include "DDSegmentation/MegatileLayerGridXY.h" diff --git a/detector/calorimeter/SEcal06_Barrel.cpp b/detector/calorimeter/SEcal06_Barrel.cpp index 3fcd03362..587ff1dc4 100644 --- a/detector/calorimeter/SEcal06_Barrel.cpp +++ b/detector/calorimeter/SEcal06_Barrel.cpp @@ -9,10 +9,8 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/SEcal06_Helpers.cpp b/detector/calorimeter/SEcal06_Helpers.cpp index a7b26afd7..10b2a7652 100644 --- a/detector/calorimeter/SEcal06_Helpers.cpp +++ b/detector/calorimeter/SEcal06_Helpers.cpp @@ -1,6 +1,5 @@ #include "SEcal06_Helpers.h" #include "LcgeoExceptions.h" -#include #include "DDSegmentation/MegatileLayerGridXY.h" #include "DDSegmentation/WaferGridXY.h" diff --git a/detector/calorimeter/SEcal06_Helpers.h b/detector/calorimeter/SEcal06_Helpers.h index 8ef75ac69..b6734dc85 100644 --- a/detector/calorimeter/SEcal06_Helpers.h +++ b/detector/calorimeter/SEcal06_Helpers.h @@ -6,7 +6,6 @@ #include "DDRec/DetectorData.h" #include "XML/Layering.h" -#include "XML/Utilities.h" #include "DD4hep/Segmentations.h" diff --git a/detector/calorimeter/SHcalSc04_Barrel_v01.cpp b/detector/calorimeter/SHcalSc04_Barrel_v01.cpp index 9ca44f8f0..51b1251e9 100644 --- a/detector/calorimeter/SHcalSc04_Barrel_v01.cpp +++ b/detector/calorimeter/SHcalSc04_Barrel_v01.cpp @@ -6,7 +6,6 @@ //==================================================================== #include "DD4hep/Printout.h" #include "DD4hep/DetFactoryHelper.h" -#include "XML/Layering.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/TiledLayerGridXY.h" diff --git a/detector/calorimeter/SHcalSc04_Barrel_v02.cpp b/detector/calorimeter/SHcalSc04_Barrel_v02.cpp index 1a82893f0..2159fcef9 100644 --- a/detector/calorimeter/SHcalSc04_Barrel_v02.cpp +++ b/detector/calorimeter/SHcalSc04_Barrel_v02.cpp @@ -7,7 +7,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" -#include "LcgeoExceptions.h" using namespace std; diff --git a/detector/calorimeter/SHcalSc04_Barrel_v03.cpp b/detector/calorimeter/SHcalSc04_Barrel_v03.cpp index 72b427b57..64ede1edf 100644 --- a/detector/calorimeter/SHcalSc04_Barrel_v03.cpp +++ b/detector/calorimeter/SHcalSc04_Barrel_v03.cpp @@ -9,9 +9,7 @@ #include "XML/Utilities.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/TiledLayerGridXY.h" -#include "LcgeoExceptions.h" -#include #include using namespace std; diff --git a/detector/calorimeter/SHcalSc04_EndcapRing.cpp b/detector/calorimeter/SHcalSc04_EndcapRing.cpp index dbf869a72..82b945de4 100644 --- a/detector/calorimeter/SHcalSc04_EndcapRing.cpp +++ b/detector/calorimeter/SHcalSc04_EndcapRing.cpp @@ -35,7 +35,6 @@ */ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "DD4hep/Shapes.h" #include "XML/Utilities.h" diff --git a/detector/calorimeter/SHcalSc04_EndcapRing_v01.cpp b/detector/calorimeter/SHcalSc04_EndcapRing_v01.cpp index c11a7b2f5..1938c06d7 100644 --- a/detector/calorimeter/SHcalSc04_EndcapRing_v01.cpp +++ b/detector/calorimeter/SHcalSc04_EndcapRing_v01.cpp @@ -38,7 +38,6 @@ */ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XML/Layering.h" #include "DD4hep/Shapes.h" #include "XML/Utilities.h" diff --git a/detector/calorimeter/Yoke05_Barrel.cpp b/detector/calorimeter/Yoke05_Barrel.cpp index cae888afc..8fcbb4e9a 100644 --- a/detector/calorimeter/Yoke05_Barrel.cpp +++ b/detector/calorimeter/Yoke05_Barrel.cpp @@ -32,7 +32,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DetType.h" #include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/calorimeter/dual-readout/include/DRconstructor.h b/detector/calorimeter/dual-readout/include/DRconstructor.h new file mode 100644 index 000000000..631576697 --- /dev/null +++ b/detector/calorimeter/dual-readout/include/DRconstructor.h @@ -0,0 +1,71 @@ +#ifndef DRconstructor_h +#define DRconstructor_h 1 + +#include "detectorSegmentations/GridDRcaloHandle_k4geo.h" + +#include "DD4hep/DetFactoryHelper.h" + +namespace ddDRcalo { + class DRconstructor { + public: + DRconstructor(xml_det_t& x_det); + ~DRconstructor() {} + + void setExpHall(dd4hep::Assembly* experimentalHall) { fExperimentalHall = experimentalHall; } + void setDRparamBarrel(dd4hep::DDSegmentation::DRparamBarrel_k4geo* paramBarrel) { fParamBarrel = paramBarrel; } + void setDRparamEndcap(dd4hep::DDSegmentation::DRparamEndcap_k4geo* paramEndcap) { fParamEndcap = paramEndcap; } + void setDescription(dd4hep::Detector* description) { fDescription = description; } + void setDetElement(dd4hep::DetElement* drDet) { fDetElement = drDet; } + void setSipmSurf(dd4hep::OpticalSurface* sipmSurf) { fSipmSurf = sipmSurf; } + void setSensDet(dd4hep::SensitiveDetector* sensDet) { + fSensDet = sensDet; + fSegmentation = dynamic_cast( sensDet->readout().segmentation().segmentation() ); + } + + void construct(); + + private: + void initiateFibers(); + void implementTowers(xml_comp_t& x_theta, dd4hep::DDSegmentation::DRparamBase_k4geo* param, dd4hep::Volume& AssemblyBoxVol); + void placeAssembly(dd4hep::DDSegmentation::DRparamBase_k4geo* param, dd4hep::Volume& AssemblyBoxVol, + dd4hep::Volume& towerVol, dd4hep::Volume& sipmWaferVol, int towerNo, int nPhi, bool isRHS=true); + void implementFibers(xml_comp_t& x_theta, dd4hep::Volume& towerVol, dd4hep::Trap& trap, dd4hep::DDSegmentation::DRparamBase_k4geo* param); + void implementFiber(dd4hep::Volume& towerVol, dd4hep::Position pos, int col, int row, float fiberLen = 200.*dd4hep::cm); + double calculateDistAtZ(TGeoTrap* rootTrap, dd4hep::Position& pos, double* norm, double z); + float calculateFiberLen(TGeoTrap* rootTrap, dd4hep::Position& pos, double* norm, double z1, double diff, double towerHeight); + dd4hep::Box calculateFullBox(TGeoTrap* rootTrap, int& rmin, int& rmax, int& cmin, int& cmax, double dz); + bool checkContained(TGeoTrap* rootTrap, dd4hep::Position& pos, double z, bool throwExcept=false); + void getNormals(TGeoTrap* rootTrap, int numxBl2, double z, double* norm1, double* norm2, double* norm3, double* norm4); + void placeUnitBox(dd4hep::Volume& fullBox, dd4hep::Volume& unitBox, int rmin, int rmax, int cmin, int cmax, bool& isEvenRow, bool& isEvenCol); + + xml_det_t fX_det; + xml_comp_t fX_barrel; + xml_comp_t fX_endcap; + xml_comp_t fX_sipmDim; + xml_comp_t fX_struct; + xml_comp_t fX_dim; + xml_comp_t fX_cladC; + xml_comp_t fX_coreC; + xml_comp_t fX_coreS; + xml_comp_t fX_worldTube; + xml_comp_t fX_barrelTube; + dd4hep::Assembly* fExperimentalHall; + dd4hep::Detector* fDescription; + dd4hep::DDSegmentation::DRparamBarrel_k4geo* fParamBarrel; + dd4hep::DDSegmentation::DRparamEndcap_k4geo* fParamEndcap; + dd4hep::DetElement* fDetElement; + dd4hep::SensitiveDetector* fSensDet; + dd4hep::OpticalSurface* fSipmSurf; + dd4hep::DDSegmentation::GridDRcalo_k4geo* fSegmentation; + + bool fVis; + int fNumx, fNumy; + std::vector< std::pair > fFiberCoords; + + std::vector< dd4hep::Tube > fFiberEnvVec; + std::vector< dd4hep::Tube > fFiberCoreCVec; + std::vector< dd4hep::Tube > fFiberCoreSVec; + }; +} + +#endif diff --git a/detector/calorimeter/dual-readout/src/DRconstructor.cpp b/detector/calorimeter/dual-readout/src/DRconstructor.cpp new file mode 100644 index 000000000..c01afa77f --- /dev/null +++ b/detector/calorimeter/dual-readout/src/DRconstructor.cpp @@ -0,0 +1,382 @@ +#include "DRconstructor.h" + +ddDRcalo::DRconstructor::DRconstructor(xml_det_t& x_det) +: fX_det(x_det), + // no default initializer for xml_comp_t + fX_barrel( x_det.child( _Unicode(barrel) ) ), + fX_endcap( x_det.child( _Unicode(endcap) ) ), + fX_sipmDim( x_det.child( _Unicode(sipmDim) ) ), + fX_struct( x_det.child( _Unicode(structure) ) ), + fX_dim( fX_struct.child( _Unicode(dim) ) ), + fX_cladC( fX_struct.child( _Unicode(cladC) ) ), + fX_coreC( fX_struct.child( _Unicode(coreC) ) ), + fX_coreS( fX_struct.child( _Unicode(coreS) ) ), + fX_worldTube( fX_struct.child( _Unicode(worldTube) ) ), + fX_barrelTube( fX_struct.child( _Unicode(barrelTube) ) ) { + fExperimentalHall = nullptr; + fParamBarrel = nullptr; + fDescription = nullptr; + fDetElement = nullptr; + fSensDet = nullptr; + fSipmSurf = nullptr; + fSegmentation = nullptr; + fVis = false; + fNumx = 0; + fNumy = 0; + fFiberCoords.reserve(100000); + fFiberEnvVec.reserve(4000); + fFiberCoreCVec.reserve(4000); + fFiberCoreSVec.reserve(4000); +} + +void ddDRcalo::DRconstructor::construct() { + // set vis on/off + fVis = fDescription->visAttributes(fX_det.visStr()).showDaughters(); + + // Tube to cover all +eta DRC geometry, Rmin = 252 - 1 mm (Rmin in endcap region), Rmax = 4500 + 1 mm, length = 4500 + 1mm + dd4hep::Tube WorldTube(fX_worldTube.rmin(), fX_worldTube.rmax(), fX_worldTube.height()/2. ); + // Tube to be subtracted from Large tube, this is for barrel region, will be subtracted from z = 0 to z = 2500 -1 mm + // Rmin = 252 - 1 mm mm, Rmax = 2500 - 1mm, length = 2500 -1 mm + dd4hep::Tube BarrelTube(fX_barrelTube.rmin(), fX_barrelTube.rmax(), fX_barrelTube.height()/2.); + // Create World assembly tube by subtracting BarrelTube from WorldTube + dd4hep::SubtractionSolid AssemblyTube(WorldTube, BarrelTube, dd4hep::Position( 0, 0, -(fX_worldTube.height() - fX_barrelTube.height())/2. )); + dd4hep::Volume AssemblyTubeVol("AssemblyTube",AssemblyTube, fDescription->material(fX_worldTube.materialStr()) ); + AssemblyTubeVol.setVisAttributes(*fDescription, fX_worldTube.visStr()); + dd4hep::PlacedVolume PlacedAssemblyTubeVol = fExperimentalHall->placeVolume( AssemblyTubeVol, dd4hep::Position( 0, 0, fX_worldTube.height()/2.) ); + PlacedAssemblyTubeVol.addPhysVolID("assembly", 0); // TODO : Have to update box ID encoding + + initiateFibers(); + + implementTowers(fX_barrel, fParamBarrel, AssemblyTubeVol); + implementTowers(fX_endcap, fParamEndcap, AssemblyTubeVol); + + if (fX_det.reflect()) { + auto refl_pos = dd4hep::Transform3D( dd4hep::RotationZYX(0., 0., M_PI), dd4hep::Position( 0, 0, -(fX_worldTube.height()/2.)) ); + dd4hep::PlacedVolume PlacedAssemblyTubeVol_refl = fExperimentalHall->placeVolume( AssemblyTubeVol, 1, refl_pos); + PlacedAssemblyTubeVol_refl.addPhysVolID("assembly", 1); + } +} + +void ddDRcalo::DRconstructor::initiateFibers() { + for (int idx = 1; idx <= 4000; idx++) { + double length = 0.05*dd4hep::cm * idx; // fiber with length interval of 0.5 mm + + fFiberEnvVec.push_back( dd4hep::Tube (0., fX_cladC.rmax(), length/2. ) ); + fFiberCoreCVec.push_back( dd4hep::Tube (0., fX_coreC.rmin(), length/2.) ); + fFiberCoreSVec.push_back( dd4hep::Tube (0., fX_coreS.rmin(), length/2.) ); + } +} + +void ddDRcalo::DRconstructor::implementTowers(xml_comp_t& x_theta, dd4hep::DDSegmentation::DRparamBase_k4geo* param, dd4hep::Volume& AssemblyBoxVol) { + double currentTheta = x_theta.theta(); + int towerNo = x_theta.start(); + for (xml_coll_t x_dThetaColl(x_theta,_U(deltatheta)); x_dThetaColl; ++x_dThetaColl, ++towerNo ) { + xml_comp_t x_deltaTheta = x_dThetaColl; + + // always use RHS for the reference + param->SetIsRHS(true); + param->SetDeltaTheta(x_deltaTheta.deltatheta()); + + double currentToC = currentTheta + x_deltaTheta.deltatheta()/2.; + currentTheta += x_deltaTheta.deltatheta(); + param->SetThetaOfCenter(currentToC); + param->init(); + + dd4hep::Trap tower( x_theta.height()/2., 0., 0., param->GetH1(), param->GetBl1(), param->GetTl1(), 0., + param->GetH2(), param->GetBl2(), param->GetTl2(), 0. ); + + dd4hep::Volume towerVol( "tower", tower, fDescription->material(x_theta.materialStr()) ); + towerVol.setVisAttributes(*fDescription, x_theta.visStr()); + + implementFibers(x_theta, towerVol, tower, param); + + xml_comp_t x_wafer ( fX_sipmDim.child( _Unicode(sipmWafer) ) ); + + // Photosensitive wafer + dd4hep::Trap sipmWaferBox( x_wafer.height()/2., 0., 0., param->GetH2(), param->GetBl2(), param->GetTl2(), 0., + param->GetH2(), param->GetBl2(), param->GetTl2(), 0. ); + dd4hep::Volume sipmWaferVol( "sipmWafer", sipmWaferBox, fDescription->material(x_wafer.materialStr()) ); + if (fVis) sipmWaferVol.setVisAttributes(*fDescription, x_wafer.visStr()); + dd4hep::SkinSurface(*fDescription, *fDetElement, "SiPMSurf_Tower"+std::to_string(towerNo), *fSipmSurf, sipmWaferVol); + + if (x_wafer.isSensitive()) { + sipmWaferVol.setSensitiveDetector(*fSensDet); + } + + // Remove sipmLayer, clear fFiberCoords instead of implementSipms() + fFiberCoords.clear(); + + for (int nPhi = 0; nPhi < x_theta.nphi(); nPhi++) { + placeAssembly(param,AssemblyBoxVol,towerVol,sipmWaferVol,towerNo,nPhi); + } + } + + param->filled(); + param->SetTotTowerNum( towerNo - x_theta.start() ); +} + +void ddDRcalo::DRconstructor::placeAssembly(dd4hep::DDSegmentation::DRparamBase_k4geo* param, dd4hep::Volume& AssemblyBoxVol, + dd4hep::Volume& towerVol, dd4hep::Volume& sipmWaferVol, int towerNo, int nPhi, bool isRHS) { + param->SetIsRHS(isRHS); + int towerNoLR = param->signedTowerNo(towerNo); + auto towerId64 = fSegmentation->setVolumeID( towerNoLR, nPhi ); + int towerId32 = fSegmentation->getFirst32bits(towerId64); + + dd4hep::Position towerPos = param->GetTowerPos(nPhi) + dd4hep::Position(0, 0, -(fX_worldTube.height()/2.)); + AssemblyBoxVol.placeVolume( towerVol, towerId32, dd4hep::Transform3D( param->GetRotationZYX(nPhi), towerPos ) ); + + // Remove sipmLayer + dd4hep::Position sipmPos = param->GetSipmLayerPos(nPhi) + dd4hep::Position(0, 0, -(fX_worldTube.height()/2.)); + dd4hep::PlacedVolume sipmWaferPhys = AssemblyBoxVol.placeVolume( sipmWaferVol, towerId32, dd4hep::Transform3D( param->GetRotationZYX(nPhi), sipmPos ) ); + sipmWaferPhys.addPhysVolID("eta", towerNoLR); + sipmWaferPhys.addPhysVolID("phi", nPhi); + sipmWaferPhys.addPhysVolID("module", 0); + + return; +} + +void ddDRcalo::DRconstructor::implementFibers(xml_comp_t& x_theta, dd4hep::Volume& towerVol, dd4hep::Trap& trap, dd4hep::DDSegmentation::DRparamBase_k4geo* param) { + auto rootTrap = trap.access(); + + float sipmSize = fX_dim.dx(); + float gridSize = fX_dim.distance(); + float towerHeight = x_theta.height(); + + float diff = fX_cladC.rmax(); // can be arbitrary small number + float z1 = towerHeight/2.-2*diff; // can be arbitrary number slightly smaller than towerHeight/2-diff + + fNumx = static_cast( std::floor( ( param->GetTl2()*2. - sipmSize/2. )/gridSize ) ) + 1; // in phi direction + fNumy = static_cast( std::floor( ( param->GetH2()*2. - sipmSize/2. )/gridSize ) ) + 1; // in eta direction + int numxBl2 = static_cast( std::floor( ( param->GetBl2()*2. - sipmSize/2. )/gridSize ) ) + 1; // only used for estimating normals + + // full length fibers + int rmin = 0, rmax = 0, cmin = 0, cmax = 0; + dd4hep::Box fullBox = calculateFullBox(rootTrap,rmin,rmax,cmin,cmax,rootTrap->GetDz()); + dd4hep::Volume fullBoxVol("fullBox",fullBox,fDescription->material(x_theta.materialStr())); + fullBoxVol.setVisAttributes(*fDescription, x_theta.visStr()); + + dd4hep::Box unitBox = dd4hep::Box(gridSize,gridSize,x_theta.height()/2.); + dd4hep::Volume unitBoxVol("unitBox",unitBox,fDescription->material(x_theta.materialStr())); + + if (fVis) + unitBoxVol.setVisAttributes(*fDescription, x_theta.visStr()); + + // // Remove cap (mirror or black paint in front of the fiber) + implementFiber(unitBoxVol, dd4hep::Position(-gridSize/2.,-gridSize/2.,0.), cmin, rmin ); + implementFiber(unitBoxVol, dd4hep::Position(gridSize/2.,-gridSize/2.,0.), cmin+1, rmin ); + implementFiber(unitBoxVol, dd4hep::Position(-gridSize/2.,gridSize/2.,0.), cmin, rmin+1 ); + implementFiber(unitBoxVol, dd4hep::Position(gridSize/2.,gridSize/2.,0.), cmin+1, rmin+1); + + + bool isEvenRow = false, isEvenCol = false; + placeUnitBox(fullBoxVol,unitBoxVol,rmin,rmax,cmin,cmax,isEvenRow,isEvenCol); + towerVol.placeVolume(fullBoxVol); + + // get normals to each side + double norm1[3] = {0.,0.,0.}, norm2[3] = {0.,0.,0.}, norm3[3] = {0.,0.,0.}, norm4[3] = {0.,0.,0.}; + getNormals(rootTrap,numxBl2,z1,norm1,norm2,norm3,norm4); + + for (int row = 0; row < fNumy; row++) { + for (int column = 0; column < fNumx; column++) { + auto localPosition = fSegmentation->localPosition(fNumx,fNumy,column,row); + dd4hep::Position pos = dd4hep::Position(localPosition); + + if ( row >= rmin && row <= rmax && column >= cmin && column <= cmax ) { + if ( ( !isEvenRow && row==rmax ) || ( !isEvenCol && column==cmax ) ) { + double pos_[3] = {pos.x(),pos.y(),-fullBox.z()+TGeoShape::Tolerance()}; + bool check = fullBox.access()->Contains(pos_); + + if (check) { + implementFiber(fullBoxVol, pos, column, row); + fFiberCoords.push_back( std::make_pair(column,row) ); + } + } + } else { + // outside tower + if (!checkContained(rootTrap,pos,z1)) continue; + + double* normX = nullptr; + double* normY = nullptr; + + // select two closest orthogonal sides + if (column > fNumx/2) normX = norm2; + else normX = norm4; + + if (row > fNumy/2) normY = norm3; + else normY = norm1; + + // compare and choose the shortest fiber length + float cand1 = calculateFiberLen(rootTrap, pos, normX, z1, diff, towerHeight); + float cand2 = calculateFiberLen(rootTrap, pos, normY, z1, diff, towerHeight); + float fiberLen = std::min(cand1,cand2); + + // not enough space to place fiber + if ( fiberLen < 0. ) continue; + + // trim fiber length in the case calculated length is longer than tower height + if (fiberLen > towerHeight) fiberLen = towerHeight; + float centerZ = towerHeight/2. - fiberLen/2.; + + // final check + if ( checkContained(rootTrap,pos,towerHeight/2.-fiberLen) ) { + dd4hep::Position centerPos( pos.x(),pos.y(),centerZ ); + implementFiber(towerVol, centerPos, column, row, fiberLen); + fFiberCoords.push_back( std::make_pair(column,row) ); + } + } + } + } +} + +// Remove cap (mirror or black paint in front of the fiber) +void ddDRcalo::DRconstructor::implementFiber(dd4hep::Volume& towerVol, dd4hep::Position pos, int col, int row, float fiberLen) { + // Don't implement fiber if the length required is shorter than 0.5 mm + if (fiberLen < 0.05*dd4hep::cm) + return; + + int fiberIdx = int( (float) fiberLen / (float) 0.05*dd4hep::cm ) - 1; // index of fiber in fiber vectors + // Actual length of fiber to be implemented, quantized in 0.5 mm unit + float approxFiberLen = 0.05*dd4hep::cm * (fiberIdx + 1); + // Fix Z position of fiber since the length of fiber can differ in [0, 0.5) mm + dd4hep::Position fixedPos = dd4hep::Position(pos.x(), pos.y(), pos.z() + (fiberLen - approxFiberLen)/2.); + + if ( fSegmentation->IsCerenkov(col,row) ) { //c fiber + dd4hep::Volume cladVol("cladC", fFiberEnvVec.at(fiberIdx), fDescription->material(fX_cladC.materialStr())); + towerVol.placeVolume( cladVol, fixedPos ); + + if (fVis) cladVol.setVisAttributes(*fDescription, fX_cladC.visStr()); // high CPU consumption! + + dd4hep::Volume coreVol("coreC", fFiberCoreCVec.at(fiberIdx), fDescription->material(fX_coreC.materialStr())); + if (fVis) coreVol.setVisAttributes(*fDescription, fX_coreC.visStr()); + cladVol.placeVolume( coreVol ); + + coreVol.setRegion(*fDescription, fX_det.regionStr()); + cladVol.setRegion(*fDescription, fX_det.regionStr()); + } else { // s fiber + dd4hep::Volume cladVol("cladS", fFiberEnvVec.at(fiberIdx), fDescription->material(fX_coreC.materialStr())); + towerVol.placeVolume( cladVol, fixedPos ); + if (fVis) cladVol.setVisAttributes(*fDescription, fX_coreC.visStr()); + + dd4hep::Volume coreVol("coreS", fFiberCoreSVec.at(fiberIdx), fDescription->material(fX_coreS.materialStr())); + if (fVis) coreVol.setVisAttributes(*fDescription, fX_coreS.visStr()); + cladVol.placeVolume( coreVol ); + + coreVol.setRegion(*fDescription, fX_det.regionStr()); + cladVol.setRegion(*fDescription, fX_det.regionStr()); + } +} + +double ddDRcalo::DRconstructor::calculateDistAtZ(TGeoTrap* rootTrap, dd4hep::Position& pos, double* norm, double z) { + double pos_[3] = {pos.x(),pos.y(),z}; + + return rootTrap->DistFromInside(pos_,norm); +} + +float ddDRcalo::DRconstructor::calculateFiberLen(TGeoTrap* rootTrap, dd4hep::Position& pos, double* norm, double z1, double diff, double towerHeight) { + float z2 = z1+diff; + float y1 = calculateDistAtZ(rootTrap,pos,norm,z1); + float y2 = calculateDistAtZ(rootTrap,pos,norm,z2); + float ymin = std::min(y1,y2); + + // return if the distance is smaller than fiber diameter + if ( ymin < 2.*fX_cladC.rmax() ) return -1.; + + // find the point where the fiber reaches a side of the tower + float slope = (y2-y1)/diff; + float y0 = (y1*z2-y2*z1)/diff; + float z = (fX_cladC.rmax()-y0)/slope; + float fiberLen = towerHeight/2. - z; + + return fiberLen; +} + +bool ddDRcalo::DRconstructor::checkContained(TGeoTrap* rootTrap, dd4hep::Position& pos, double z, bool throwExcept) { + double pos_[3] = {pos.x(),pos.y(),z}; + bool check = rootTrap->Contains(pos_); + + if ( throwExcept && !check ) throw std::runtime_error("Fiber must be in the tower!"); + return check; +} + +void ddDRcalo::DRconstructor::getNormals(TGeoTrap* rootTrap, int numxBl2, double z, double* norm1, double* norm2, double* norm3, double* norm4) { + dd4hep::Position pos1 = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,fNumx/2,0) ); + dd4hep::Position pos2 = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,fNumx/2+numxBl2/2-1,fNumy/2) ); + dd4hep::Position pos3 = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,fNumx/2,fNumy-1) ); + dd4hep::Position pos4 = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,fNumx/2-numxBl2/2+1,fNumy/2) ); + double pos1_[3] = {pos1.x(),pos1.y(),z}; + double pos2_[3] = {pos2.x(),pos2.y(),z}; + double pos3_[3] = {pos3.x(),pos3.y(),z}; + double pos4_[3] = {pos4.x(),pos4.y(),z}; + double dir[3] = {0.,0.,0.}; + + rootTrap->ComputeNormal(pos1_,dir,norm1); + rootTrap->ComputeNormal(pos2_,dir,norm2); + rootTrap->ComputeNormal(pos3_,dir,norm3); + rootTrap->ComputeNormal(pos4_,dir,norm4); + norm1[2] = 0.; // check horizontal distance only + norm2[2] = 0.; + norm3[2] = 0.; + norm4[2] = 0.; +} + +dd4hep::Box ddDRcalo::DRconstructor::calculateFullBox(TGeoTrap* rootTrap, int& rmin, int& rmax, int& cmin, int& cmax, double dz) { + float gridSize = fX_dim.distance(); + double zmin = -rootTrap->GetDz() + TGeoShape::Tolerance(); + float xmin = 0., xmax = 0., ymin = 0., ymax = 0.; + + for (int row = 0; row < fNumy; row++) { // bottom-up + auto localPosition = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,fNumx/2,row) ); + auto pos = localPosition + dd4hep::Position(0.,-gridSize/2.,0.); + if ( checkContained(rootTrap,pos,zmin) ) { + ymin = pos.y(); + rmin = row; + break; + } + } + + for (int row = fNumy-1; row !=0 ; row--) { // top-down + auto localPosition = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,fNumx/2,row) ); + auto pos = localPosition + dd4hep::Position(0.,gridSize/2.,0.); + if ( checkContained(rootTrap,pos,zmin) ) { + ymax = pos.y(); + rmax = row; + break; + } + } + + for (int col = 0; col < fNumx; col++) { // left-right + auto localPosition = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,col,rmin) ); + auto pos = localPosition + dd4hep::Position(-gridSize/2.,-gridSize/2.,0.); + if ( checkContained(rootTrap,pos,zmin) ) { + xmin = pos.x(); + cmin = col; + break; + } + } + + for (int col = fNumx-1; col!=0; col--) { // right-left + auto localPosition = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,col,rmin) ); + auto pos = localPosition + dd4hep::Position(gridSize/2.,-gridSize/2.,0.); + if ( checkContained(rootTrap,pos,zmin) ) { + xmax = pos.x(); + cmax = col; + break; + } + } + + return dd4hep::Box( (xmax-xmin)/2., (ymax-ymin)/2., dz ); +} + +void ddDRcalo::DRconstructor::placeUnitBox(dd4hep::Volume& fullBox, dd4hep::Volume& unitBox, int rmin, int rmax, int cmin, int cmax, bool& isEvenRow, bool& isEvenCol) { + for (int row = rmin; row < rmax; row+=2) { + for (int col = cmin; col < cmax; col+=2) { + auto pos0 = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,col,row) ); + auto pos3 = dd4hep::Position( fSegmentation->localPosition(fNumx,fNumy,col+1,row+1) ); + fullBox.placeVolume(unitBox,(pos0+pos3)/2.); + } + } + + isEvenRow = (rmax-rmin+1)%2==0; + isEvenCol = (cmax-cmin+1)%2==0; + return; +} diff --git a/detector/calorimeter/dual-readout/src/FiberDualReadoutCalo_o1_v01.cpp b/detector/calorimeter/dual-readout/src/FiberDualReadoutCalo_o1_v01.cpp new file mode 100644 index 000000000..11140c0a5 --- /dev/null +++ b/detector/calorimeter/dual-readout/src/FiberDualReadoutCalo_o1_v01.cpp @@ -0,0 +1,72 @@ +#include "detectorSegmentations/DRparamBarrel_k4geo.h" +#include "detectorSegmentations/DRparamEndcap_k4geo.h" + +#include "DRconstructor.h" + +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/OpticalSurfaces.h" +#include "DD4hep/Detector.h" + +namespace ddDRcalo { + static dd4hep::Ref_t create_detector( dd4hep::Detector &description, xml_h xmlElement, dd4hep::SensitiveDetector sensDet ) { + + // Get the detector description from the xml-tree + xml_det_t x_det = xmlElement; + std::string name = x_det.nameStr(); + // Create the detector element + dd4hep::DetElement drDet( name, x_det.id() ); + // set the sensitive detector type to the DD4hep calorimeter + dd4hep::xml::Dimension sensDetType = xmlElement.child(_Unicode(sensitive)); + sensDet.setType(sensDetType.typeStr()); + // Get the world volume + dd4hep::Assembly experimentalHall("hall"); + // Get the dimensions defined in the xml-tree + xml_comp_t x_barrel ( x_det.child( _Unicode(barrel) ) ); + xml_comp_t x_endcap ( x_det.child( _Unicode(endcap) ) ); + xml_comp_t x_structure ( x_det.child( _Unicode(structure) ) ); + xml_comp_t x_dim ( x_structure.child( _Unicode(dim) ) ); + xml_comp_t x_sipmDim ( x_det.child( _Unicode(sipmDim) ) ); + + dd4hep::OpticalSurfaceManager surfMgr = description.surfaceManager(); + dd4hep::OpticalSurface sipmSurfProp = surfMgr.opticalSurface("/world/"+name+"#SiPMSurf"); + surfMgr.opticalSurface("/world/"+name+"#FilterSurf"); // actual filtering applied in the stepping action + + auto segmentation = dynamic_cast( sensDet.readout().segmentation().segmentation() ); + segmentation->setGridSize( x_dim.distance() ); + segmentation->setSipmSize( x_dim.dx() ); + + auto paramBarrel = segmentation->paramBarrel(); + paramBarrel->SetInnerX(x_barrel.rmin()); + paramBarrel->SetTowerH(x_barrel.height()); + paramBarrel->SetNumZRot(x_barrel.nphi()); + paramBarrel->SetSipmHeight(x_sipmDim.height()); + + auto paramEndcap = segmentation->paramEndcap(); + paramEndcap->SetInnerX(x_endcap.rmin()); + paramEndcap->SetTowerH(x_endcap.height()); + paramEndcap->SetNumZRot(x_endcap.nphi()); + paramEndcap->SetSipmHeight(x_sipmDim.height()); + + auto constructor = DRconstructor(x_det); + constructor.setExpHall(&experimentalHall); + constructor.setDRparamBarrel(paramBarrel); + constructor.setDRparamEndcap(paramEndcap); + constructor.setDescription(&description); + constructor.setDetElement(&drDet); + constructor.setSipmSurf(&sipmSurfProp); + constructor.setSensDet(&sensDet); + constructor.construct(); // right + + dd4hep::Volume worldVol = description.pickMotherVolume(drDet); + dd4hep::PlacedVolume hallPlace = worldVol.placeVolume(experimentalHall); + hallPlace.addPhysVolID("system",x_det.id()); + // connect placed volume and physical volume + drDet.setPlacement( hallPlace ); + + paramBarrel->finalized(); + paramEndcap->finalized(); + + return drDet; + } +} // namespace detector +DECLARE_DETELEMENT(FiberDualReadoutCalo_o1_v01, ddDRcalo::create_detector) // factory method diff --git a/detector/fcal/BeamCal_o1_v01_geo.cpp b/detector/fcal/BeamCal_o1_v01_geo.cpp index c5f660f5e..b7bef1127 100644 --- a/detector/fcal/BeamCal_o1_v01_geo.cpp +++ b/detector/fcal/BeamCal_o1_v01_geo.cpp @@ -1,5 +1,4 @@ #include -#include "DD4hep/DetType.h" #include #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/fcal/BeamCal_o1_v02_geo.cpp b/detector/fcal/BeamCal_o1_v02_geo.cpp index 3721299b7..53196e85f 100644 --- a/detector/fcal/BeamCal_o1_v02_geo.cpp +++ b/detector/fcal/BeamCal_o1_v02_geo.cpp @@ -1,5 +1,4 @@ #include -#include "DD4hep/DetType.h" #include #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/muonSystem/muonSystemMuRWELL_o1_v01.cpp b/detector/muonSystem/muonSystemMuRWELL_o1_v01.cpp index 539f6f8e6..2f4321433 100644 --- a/detector/muonSystem/muonSystemMuRWELL_o1_v01.cpp +++ b/detector/muonSystem/muonSystemMuRWELL_o1_v01.cpp @@ -23,10 +23,8 @@ If used with sensitive layers, the readout must contain a "slice" field #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "XML/XMLElements.h" -#include -#include +#include #include using namespace std; @@ -526,7 +524,7 @@ static dd4hep::Ref_t createmuonSystemMuRWELL_o1_v01(dd4hep::Detector& lcdd, dd4hep::Position detectorLayerTrans(0., 0., 0.); dd4hep::PlacedVolume detectorLayerPhys = BarrelVolume.placeVolume(BarrelDetectorLayerVolume, dd4hep::Transform3D(dd4hep::RotationZ(0.), detectorLayerTrans)); detectorLayerPhys.addPhysVolID("layer", numBarrelLayer+1); - dd4hep::DetElement BarrelDetectorLayerDE(BarrelDE, "MSBarrelDetectorLayerDE" + numBarrelLayer+1, numBarrelLayer+1); + dd4hep::DetElement BarrelDetectorLayerDE(BarrelDE, "MSBarrel_DetectorLayerDE_" + std::to_string(numBarrelLayer+1), numBarrelLayer+1); BarrelDetectorLayerDE.setPlacement(detectorLayerPhys); BarrelDetectorLayerVolume.setVisAttributes(lcdd.visAttributes("no_vis")); @@ -591,7 +589,7 @@ static dd4hep::Ref_t createmuonSystemMuRWELL_o1_v01(dd4hep::Detector& lcdd, dd4hep::Position radiatorLayerTrans(0., 0., 0.); dd4hep::PlacedVolume radiatorLayerPhys = BarrelVolume.placeVolume(BarrelRadiatorLayerVolume, dd4hep::Transform3D(dd4hep::RotationZ(0.), radiatorLayerTrans)); - dd4hep::DetElement BarrelRadiatorLayerDE(BarrelDE, "MSBarrel_RadiatorLayerDE" + numBarrelRadiatorLayer+1, numBarrelRadiatorLayer+1); + dd4hep::DetElement BarrelRadiatorLayerDE(BarrelDE, "MSBarrel_RadiatorLayerDE_" + std::to_string(numBarrelRadiatorLayer+1), numBarrelRadiatorLayer+1); BarrelRadiatorLayerDE.setPlacement(radiatorLayerPhys); BarrelRadiatorLayerVolume.setVisAttributes(lcdd.visAttributes("no_vis")); @@ -616,7 +614,7 @@ static dd4hep::Ref_t createmuonSystemMuRWELL_o1_v01(dd4hep::Detector& lcdd, endcapTrans = dd4hep::Position(0., 0., endcapType * endcapOffset); endcapPhys = detectorVolume.placeVolume(endcapVolume, dd4hep::Transform3D(dd4hep::RotationZ(0.), endcapTrans)); endcapPhys.addPhysVolID("type", endcapType); - EndcapDE = dd4hep::DetElement(detElement, "MSEndcapDE" + endcapType , endcapType); + EndcapDE = dd4hep::DetElement(detElement, "MSEndcapDE_" + std::to_string(endcapType) , endcapType); EndcapDE.setPlacement(endcapPhys); endcapVolume.setVisAttributes(lcdd.visAttributes("no_vis")); diff --git a/detector/other/Beampipe_o1_v01_geo.cpp b/detector/other/Beampipe_o1_v01_geo.cpp index 676275a83..b30f77b30 100644 --- a/detector/other/Beampipe_o1_v01_geo.cpp +++ b/detector/other/Beampipe_o1_v01_geo.cpp @@ -9,14 +9,10 @@ #include "OtherDetectorHelpers.h" #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" #include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include diff --git a/detector/other/BoxSupport_o1_v01_geo.cpp b/detector/other/BoxSupport_o1_v01_geo.cpp index 1c3ec3e54..204c0d3a2 100644 --- a/detector/other/BoxSupport_o1_v01_geo.cpp +++ b/detector/other/BoxSupport_o1_v01_geo.cpp @@ -1,11 +1,7 @@ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include using namespace std; diff --git a/detector/other/ConicalSupport_o1_v01_geo.cpp b/detector/other/ConicalSupport_o1_v01_geo.cpp index a5462e394..7dd25b937 100644 --- a/detector/other/ConicalSupport_o1_v01_geo.cpp +++ b/detector/other/ConicalSupport_o1_v01_geo.cpp @@ -1,11 +1,8 @@ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include using namespace std; diff --git a/detector/other/FieldMapBrBz.cpp b/detector/other/FieldMapBrBz.cpp index 1a6036252..e14b92519 100644 --- a/detector/other/FieldMapBrBz.cpp +++ b/detector/other/FieldMapBrBz.cpp @@ -7,11 +7,9 @@ #include #endif -#include #include -#include #include #include diff --git a/detector/other/FieldMapXYZ.cpp b/detector/other/FieldMapXYZ.cpp index 9da75be12..8857cff18 100644 --- a/detector/other/FieldMapXYZ.cpp +++ b/detector/other/FieldMapXYZ.cpp @@ -7,14 +7,10 @@ #include #endif -#include #include -#include -//Thanks stackoverflow -#include #include diff --git a/detector/other/Mask_o1_v01_geo.cpp b/detector/other/Mask_o1_v01_geo.cpp index 16ebefdfb..08df477ad 100644 --- a/detector/other/Mask_o1_v01_geo.cpp +++ b/detector/other/Mask_o1_v01_geo.cpp @@ -10,7 +10,6 @@ #include "OtherDetectorHelpers.h" #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include #include diff --git a/detector/other/MaterialEnvelope_o1_v01.cpp b/detector/other/MaterialEnvelope_o1_v01.cpp index e3e26b878..8513803f3 100644 --- a/detector/other/MaterialEnvelope_o1_v01.cpp +++ b/detector/other/MaterialEnvelope_o1_v01.cpp @@ -7,7 +7,6 @@ // //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "XML/Layering.h" #include "XML/Utilities.h" using namespace std; diff --git a/detector/other/PolyhedralBarrelSurfaces_geo.cpp b/detector/other/PolyhedralBarrelSurfaces_geo.cpp index 72d9c6b5f..cf270f88e 100644 --- a/detector/other/PolyhedralBarrelSurfaces_geo.cpp +++ b/detector/other/PolyhedralBarrelSurfaces_geo.cpp @@ -11,7 +11,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "DDRec/Surface.h" -#include "DDRec/DetectorData.h" #include "XML/Utilities.h" using dd4hep::BUILD_ENVELOPE; diff --git a/detector/other/PolyhedralEndcapSurfaces_geo.cpp b/detector/other/PolyhedralEndcapSurfaces_geo.cpp index 63fbe8fdf..a4047ecca 100644 --- a/detector/other/PolyhedralEndcapSurfaces_geo.cpp +++ b/detector/other/PolyhedralEndcapSurfaces_geo.cpp @@ -11,7 +11,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "DDRec/Surface.h" -#include "DDRec/DetectorData.h" #include "XML/Utilities.h" using dd4hep::BUILD_ENVELOPE; diff --git a/detector/other/SCoil02_geo.cpp b/detector/other/SCoil02_geo.cpp index 970818eee..1d6a8f23c 100644 --- a/detector/other/SCoil02_geo.cpp +++ b/detector/other/SCoil02_geo.cpp @@ -5,12 +5,7 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" -#include "DDRec/Surface.h" -#include "XMLHandlerDB.h" #include "XML/Utilities.h" -#include #include "DDRec/DetectorData.h" //#include "GearWrapper.h" diff --git a/detector/other/SServices00.cpp b/detector/other/SServices00.cpp index 6c8e2babc..ffacf584d 100644 --- a/detector/other/SServices00.cpp +++ b/detector/other/SServices00.cpp @@ -28,9 +28,7 @@ // #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include "XML/Utilities.h" -#include "DD4hep/DetType.h" #include "DDRec/Surface.h" #include "XMLHandlerDB.h" diff --git a/detector/other/SServices00.h b/detector/other/SServices00.h index c9637c00c..7a2eb5945 100644 --- a/detector/other/SServices00.h +++ b/detector/other/SServices00.h @@ -2,7 +2,6 @@ #include #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include "DD4hep/DetElement.h" //===================================== diff --git a/detector/other/Solenoid_o1_v01_gep.cpp b/detector/other/Solenoid_o1_v01_gep.cpp index 89408c442..e8cb851d0 100644 --- a/detector/other/Solenoid_o1_v01_gep.cpp +++ b/detector/other/Solenoid_o1_v01_gep.cpp @@ -1,7 +1,4 @@ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" -#include "XML/Layering.h" -#include "TGeoTrd2.h" #include "XML/Utilities.h" #include "DDRec/DetectorData.h" diff --git a/detector/other/TrackerBarrelSupport_o1_v01.cpp b/detector/other/TrackerBarrelSupport_o1_v01.cpp index 3d562d040..5fe694fe6 100644 --- a/detector/other/TrackerBarrelSupport_o1_v01.cpp +++ b/detector/other/TrackerBarrelSupport_o1_v01.cpp @@ -1,13 +1,9 @@ -#include "OtherDetectorHelpers.h" #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include using namespace std; diff --git a/detector/other/TrackerEndcapSupport_o1_v01.cpp b/detector/other/TrackerEndcapSupport_o1_v01.cpp index bbd4014e2..1bc9b680c 100644 --- a/detector/other/TrackerEndcapSupport_o1_v01.cpp +++ b/detector/other/TrackerEndcapSupport_o1_v01.cpp @@ -1,13 +1,8 @@ -#include "OtherDetectorHelpers.h" #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include using namespace std; diff --git a/detector/other/TrackerEndcapSupport_o1_v02.cpp b/detector/other/TrackerEndcapSupport_o1_v02.cpp index 1f2d4d746..8d3e58111 100644 --- a/detector/other/TrackerEndcapSupport_o1_v02.cpp +++ b/detector/other/TrackerEndcapSupport_o1_v02.cpp @@ -4,16 +4,11 @@ // * an issue that reflected volumes had double thickness // * an issue that the whole layer volume was twice as thick as it needed to be -#include "OtherDetectorHelpers.h" #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include using namespace std; diff --git a/detector/other/TubeSupport_o1_v01_geo.cpp b/detector/other/TubeSupport_o1_v01_geo.cpp index cd463688b..42c916550 100644 --- a/detector/other/TubeSupport_o1_v01_geo.cpp +++ b/detector/other/TubeSupport_o1_v01_geo.cpp @@ -1,11 +1,7 @@ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DDRec/DetectorData.h" #include "DDRec/Surface.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include -#include #include using namespace std; diff --git a/detector/other/TubeX01_geo.cpp b/detector/other/TubeX01_geo.cpp index 04d0d11e0..7546eb6f3 100644 --- a/detector/other/TubeX01_geo.cpp +++ b/detector/other/TubeX01_geo.cpp @@ -5,7 +5,6 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include "DDRec/DetectorData.h" #include "XMLHandlerDB.h" #include diff --git a/detector/tracker/DriftChamber_o1_v01.cpp b/detector/tracker/DriftChamber_o1_v01.cpp index 6023596f6..23f3cb184 100644 --- a/detector/tracker/DriftChamber_o1_v01.cpp +++ b/detector/tracker/DriftChamber_o1_v01.cpp @@ -4,13 +4,8 @@ \***************************************************************************************/ #include "DD4hep/DetFactoryHelper.h" #include "DD4hep/Printout.h" -#include "DD4hep/detail/DetectorInterna.h" -#include "TClass.h" #include "TMath.h" #include "XML/Utilities.h" -#include -#include -#include using namespace std; using namespace dd4hep; diff --git a/detector/tracker/DriftChamber_o1_v02.cpp b/detector/tracker/DriftChamber_o1_v02.cpp index 87c54e0fd..a174f804f 100644 --- a/detector/tracker/DriftChamber_o1_v02.cpp +++ b/detector/tracker/DriftChamber_o1_v02.cpp @@ -15,7 +15,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "DD4hep/Shapes.h" -#include "DD4hep/Printout.h" #include "DD4hep/Detector.h" #include "DDRec/DCH_info.h" diff --git a/detector/tracker/FTD_Simple_Staggered_geo.cpp b/detector/tracker/FTD_Simple_Staggered_geo.cpp index 52049c85e..d745f73dd 100644 --- a/detector/tracker/FTD_Simple_Staggered_geo.cpp +++ b/detector/tracker/FTD_Simple_Staggered_geo.cpp @@ -5,9 +5,7 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DetType.h" #include "XMLHandlerDB.h" #include "XML/Utilities.h" diff --git a/detector/tracker/SET_Simple_Planar_geo.cpp b/detector/tracker/SET_Simple_Planar_geo.cpp index 91e906c81..c295104a3 100644 --- a/detector/tracker/SET_Simple_Planar_geo.cpp +++ b/detector/tracker/SET_Simple_Planar_geo.cpp @@ -5,8 +5,6 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" #include "DDRec/Surface.h" #include "DDRec/DetectorData.h" #include "XMLHandlerDB.h" diff --git a/detector/tracker/SIT_Simple_Pixel_geo.cpp b/detector/tracker/SIT_Simple_Pixel_geo.cpp index 01bd3572b..dc26db6ab 100644 --- a/detector/tracker/SIT_Simple_Pixel_geo.cpp +++ b/detector/tracker/SIT_Simple_Pixel_geo.cpp @@ -5,8 +5,6 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" #include "DDRec/Surface.h" #include "DDRec/DetectorData.h" #include "XML/Utilities.h" diff --git a/detector/tracker/SIT_Simple_Planar_geo.cpp b/detector/tracker/SIT_Simple_Planar_geo.cpp index e6a93e80c..71475c27f 100644 --- a/detector/tracker/SIT_Simple_Planar_geo.cpp +++ b/detector/tracker/SIT_Simple_Planar_geo.cpp @@ -5,8 +5,6 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" #include "DDRec/Surface.h" #include "DDRec/DetectorData.h" #include "XML/Utilities.h" diff --git a/detector/tracker/SiTrackerEndcap_o2_v02_geo.cpp b/detector/tracker/SiTrackerEndcap_o2_v02_geo.cpp index 034f276a0..b82e21e61 100644 --- a/detector/tracker/SiTrackerEndcap_o2_v02_geo.cpp +++ b/detector/tracker/SiTrackerEndcap_o2_v02_geo.cpp @@ -18,7 +18,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Utilities.h" #include -#include "DDRec/DetectorData.h" using namespace std; using dd4hep::Ref_t; diff --git a/detector/tracker/SiTrackerEndcap_o2_v02ext_geo.cpp b/detector/tracker/SiTrackerEndcap_o2_v02ext_geo.cpp index 37b26cac6..f96ad7d5e 100644 --- a/detector/tracker/SiTrackerEndcap_o2_v02ext_geo.cpp +++ b/detector/tracker/SiTrackerEndcap_o2_v02ext_geo.cpp @@ -21,9 +21,7 @@ #include "XML/Utilities.h" #include #include "DDRec/DetectorData.h" -#include "XML/DocumentHandler.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/SiTrackerEndcap_o2_v03_geo.cpp b/detector/tracker/SiTrackerEndcap_o2_v03_geo.cpp index 0aa30554e..0f8c221e3 100644 --- a/detector/tracker/SiTrackerEndcap_o2_v03_geo.cpp +++ b/detector/tracker/SiTrackerEndcap_o2_v03_geo.cpp @@ -18,7 +18,6 @@ #include "DD4hep/DetFactoryHelper.h" #include "XML/Utilities.h" #include -#include "DDRec/DetectorData.h" using namespace std; diff --git a/detector/tracker/TPC10_geo.cpp b/detector/tracker/TPC10_geo.cpp index 888aff300..57c911436 100644 --- a/detector/tracker/TPC10_geo.cpp +++ b/detector/tracker/TPC10_geo.cpp @@ -5,8 +5,6 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" #include "LcgeoExceptions.h" #include "lcgeo.h" #include "DDRec/Surface.h" diff --git a/detector/tracker/TrackerBarrel_o1_v03_geo.cpp b/detector/tracker/TrackerBarrel_o1_v03_geo.cpp index 07438043e..8e0b50403 100644 --- a/detector/tracker/TrackerBarrel_o1_v03_geo.cpp +++ b/detector/tracker/TrackerBarrel_o1_v03_geo.cpp @@ -14,7 +14,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/TrackerBarrel_o1_v04_geo.cpp b/detector/tracker/TrackerBarrel_o1_v04_geo.cpp index bf15cab2a..11b660450 100644 --- a/detector/tracker/TrackerBarrel_o1_v04_geo.cpp +++ b/detector/tracker/TrackerBarrel_o1_v04_geo.cpp @@ -14,7 +14,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/TrackerBarrel_o1_v05_geo.cpp b/detector/tracker/TrackerBarrel_o1_v05_geo.cpp index 317e03331..ea00cdd24 100644 --- a/detector/tracker/TrackerBarrel_o1_v05_geo.cpp +++ b/detector/tracker/TrackerBarrel_o1_v05_geo.cpp @@ -17,7 +17,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/TrackerEndcap_o1_v05_geo.cpp b/detector/tracker/TrackerEndcap_o1_v05_geo.cpp index 6e3b892b6..02ca1d324 100644 --- a/detector/tracker/TrackerEndcap_o1_v05_geo.cpp +++ b/detector/tracker/TrackerEndcap_o1_v05_geo.cpp @@ -26,7 +26,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/TrackerEndcap_o2_v05_geo.cpp b/detector/tracker/TrackerEndcap_o2_v05_geo.cpp index fa12598e5..f121df6ba 100644 --- a/detector/tracker/TrackerEndcap_o2_v05_geo.cpp +++ b/detector/tracker/TrackerEndcap_o2_v05_geo.cpp @@ -26,7 +26,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/TrackerEndcap_o2_v06_geo.cpp b/detector/tracker/TrackerEndcap_o2_v06_geo.cpp index 49d65d581..14ce44f68 100644 --- a/detector/tracker/TrackerEndcap_o2_v06_geo.cpp +++ b/detector/tracker/TrackerEndcap_o2_v06_geo.cpp @@ -30,7 +30,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include @@ -184,11 +183,11 @@ static Ref_t create_detector(Detector& theDetector, xml_h e, SensitiveDetector s Box mod_shape(m_vol.solid()); - if(r - mod_shape->GetDZ()GetDZ(); + if(r < innerR) + innerR = r; - if(r + mod_shape->GetDZ()>outerR) - outerR = r + mod_shape->GetDZ(); + if(r + 2*mod_shape->GetDY() > outerR) + outerR = r + 2*mod_shape->GetDY(); sumZ+=zstart; diff --git a/detector/tracker/VXD04_geo.cpp b/detector/tracker/VXD04_geo.cpp index 9b13b09a4..11cf56a24 100644 --- a/detector/tracker/VXD04_geo.cpp +++ b/detector/tracker/VXD04_geo.cpp @@ -5,8 +5,6 @@ // $Id$ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" #include "DDRec/Surface.h" #include "DDRec/DetectorData.h" diff --git a/detector/tracker/VertexBarrel_detailed_o1_v01_geo.cpp b/detector/tracker/VertexBarrel_detailed_o1_v01_geo.cpp index cc9704eb5..801c854c2 100644 --- a/detector/tracker/VertexBarrel_detailed_o1_v01_geo.cpp +++ b/detector/tracker/VertexBarrel_detailed_o1_v01_geo.cpp @@ -3,8 +3,8 @@ // Based on ZPlanarTracker module from F. Gaede // Tracking detector to describe the FCC-ee IDEA vertex detector barrel. -// The vertex detector is assembled of stave structures which feature -// support and readout (flex) elements. Each stave features multiple +// The vertex detector is assembled of stave structures which feature +// support and readout (flex) elements. Each stave features multiple // individual modules, that consist of sensitive and insensitive // sensor elements. //-------------------------------------------------------------------- @@ -14,16 +14,9 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" #include "XML/Utilities.h" -#include "XMLHandlerDB.h" #include "DDRec/Surface.h" -#include "DDRec/DetectorData.h" -#include -#include -#include -#include "UTIL/LCTrackerConf.h" -#include using namespace std; diff --git a/detector/tracker/VertexEndcap_detailed_o1_v01_geo.cpp b/detector/tracker/VertexEndcap_detailed_o1_v01_geo.cpp index 01dc84a03..8d423703f 100644 --- a/detector/tracker/VertexEndcap_detailed_o1_v01_geo.cpp +++ b/detector/tracker/VertexEndcap_detailed_o1_v01_geo.cpp @@ -17,10 +17,8 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" -#include +#include #include "XML/Utilities.h" -#include "DD4hep/Printout.h" -#include "DDRec/DetectorData.h" using namespace std; diff --git a/detector/tracker/VertexEndcap_o1_v05_geo.cpp b/detector/tracker/VertexEndcap_o1_v05_geo.cpp index 3cd66e2ca..bc4018218 100644 --- a/detector/tracker/VertexEndcap_o1_v05_geo.cpp +++ b/detector/tracker/VertexEndcap_o1_v05_geo.cpp @@ -20,7 +20,6 @@ #include "DDRec/DetectorData.h" #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/ZPlanarTracker_geo.cpp b/detector/tracker/ZPlanarTracker_geo.cpp index bb7b5e2a0..174237a77 100644 --- a/detector/tracker/ZPlanarTracker_geo.cpp +++ b/detector/tracker/ZPlanarTracker_geo.cpp @@ -13,10 +13,8 @@ #include "DDRec/Surface.h" #include "DDRec/DetectorData.h" -#include #include -#include #include "UTIL/LCTrackerConf.h" #include diff --git a/detector/tracker/parametrised_DriftChamber_o1_v01.cpp b/detector/tracker/parametrised_DriftChamber_o1_v01.cpp index 16e28a490..969424c8f 100644 --- a/detector/tracker/parametrised_DriftChamber_o1_v01.cpp +++ b/detector/tracker/parametrised_DriftChamber_o1_v01.cpp @@ -1,8 +1,6 @@ #include "DD4hep/DetFactoryHelper.h" -#include "DD4hep/Printout.h" #include "DDSegmentation/Segmentation.h" #include "detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h" -#include "XML/Utilities.h" #include "XML/XMLElements.h" namespace det { diff --git a/detectorCommon/CMakeLists.txt b/detectorCommon/CMakeLists.txt index 0982e0a4c..407027b93 100644 --- a/detectorCommon/CMakeLists.txt +++ b/detectorCommon/CMakeLists.txt @@ -4,6 +4,7 @@ file(GLOB sources src/*.cpp) add_dd4hep_plugin(detectorCommon SHARED ${sources}) +add_library(k4geo::detectorCommon ALIAS detectorCommon) target_link_libraries(detectorCommon DD4hep::DDCore DD4hep::DDG4 detectorSegmentations) target_include_directories(detectorCommon PUBLIC diff --git a/detectorCommon/include/detectorCommon/DetUtils_k4geo.h b/detectorCommon/include/detectorCommon/DetUtils_k4geo.h index 04c1f93e6..9085a5820 100644 --- a/detectorCommon/include/detectorCommon/DetUtils_k4geo.h +++ b/detectorCommon/include/detectorCommon/DetUtils_k4geo.h @@ -16,13 +16,11 @@ #include "DDSegmentation/CartesianGridXYZ.h" #include "DDSegmentation/PolarGridRPhi.h" -// Geant #include "G4Step.hh" // CLHEP #include "CLHEP/Vector/ThreeVector.h" -#include "TGeoManager.h" /** Given a XML element with several daughters with the same name, e.g. diff --git a/detectorCommon/src/xtalk_neighbors_moduleThetaMergedSegmentation.cpp b/detectorCommon/src/xtalk_neighbors_moduleThetaMergedSegmentation.cpp index 7e4ce2397..cc065cac5 100644 --- a/detectorCommon/src/xtalk_neighbors_moduleThetaMergedSegmentation.cpp +++ b/detectorCommon/src/xtalk_neighbors_moduleThetaMergedSegmentation.cpp @@ -1,13 +1,8 @@ -#include "detectorCommon/DetUtils_k4geo.h" #include "detectorCommon/xtalk_neighbors_moduleThetaMergedSegmentation.h" -// DD4hep -#include "DDG4/Geant4Mapping.h" -#include "DDG4/Geant4VolumeManager.h" #include -#include namespace det { namespace crosstalk { diff --git a/detectorSegmentations/CMakeLists.txt b/detectorSegmentations/CMakeLists.txt index c41cbb8f9..0d7a5fb3a 100644 --- a/detectorSegmentations/CMakeLists.txt +++ b/detectorSegmentations/CMakeLists.txt @@ -5,6 +5,7 @@ file(GLOB sources src/*.cpp) add_library(detectorSegmentations SHARED ${sources} ) +add_library(k4geo::detectorSegmentations ALIAS detectorSegmentations) dd4hep_generate_rootmap(detectorSegmentationsPlugin) target_link_libraries(detectorSegmentations DD4hep::DDCore) target_include_directories(detectorSegmentations diff --git a/detectorSegmentations/include/detectorSegmentations/DRparamBarrel_k4geo.h b/detectorSegmentations/include/detectorSegmentations/DRparamBarrel_k4geo.h new file mode 100644 index 000000000..6958d6d14 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/DRparamBarrel_k4geo.h @@ -0,0 +1,24 @@ +#ifndef DETSEGMENTATION_DRPARAMBARREL_H +#define DETSEGMENTATION_DRPARAMBARREL_H + +#include "detectorSegmentations/DRparamBase_k4geo.h" + +#include "DD4hep/DetFactoryHelper.h" + + +namespace dd4hep { +namespace DDSegmentation { + class DRparamBarrel_k4geo : public DRparamBase_k4geo { + public: + DRparamBarrel_k4geo(); + virtual ~DRparamBarrel_k4geo(); + + virtual void SetDeltaThetaByTowerNo(int signedTowerNo, int) override; + virtual void SetThetaOfCenterByTowerNo(int signedTowerNo, int) override; + + virtual void init() override; + }; +} +} + +#endif diff --git a/detectorSegmentations/include/detectorSegmentations/DRparamBase_k4geo.h b/detectorSegmentations/include/detectorSegmentations/DRparamBase_k4geo.h new file mode 100644 index 000000000..c072eabb5 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/DRparamBase_k4geo.h @@ -0,0 +1,99 @@ +#ifndef DETSEGMENTATION_DRPARAMBASE_H +#define DETSEGMENTATION_DRPARAMBASE_H + +#include "TVector3.h" +#include "DD4hep/DetFactoryHelper.h" + +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + class DRparamBase_k4geo { + public: + DRparamBase_k4geo(); + virtual ~DRparamBase_k4geo(); + + void SetIsRHS(bool isRHS) { fIsRHS = isRHS; } + void SetInnerX(double innerX) { fInnerX = innerX; } + void SetTowerH(double towerH) { fTowerH = towerH; } + void SetNumZRot(int num) { fNumZRot = num; fPhiZRot = 2*M_PI/(double)num; } + void SetDeltaTheta(double theta) { fDeltaTheta = theta; } + void SetThetaOfCenter(double theta) { fThetaOfCenter = theta; } + void SetSipmHeight(double SipmHeight) { fSipmHeight = SipmHeight; } + + bool GetIsRHS() { return fIsRHS; } + double GetCurrentInnerR() { return fCurrentInnerR; } + double GetTowerH() { return fTowerH; } + double GetSipmHeight() { return fSipmHeight; } + double GetH1() { return fCurrentInnerHalf; } + double GetBl1() { return fV3.X()*std::tan(fPhiZRot/2.); } + double GetTl1() { return fV1.X()*std::tan(fPhiZRot/2.); } + double GetH2() { return fCurrentOuterHalf; } + double GetBl2() { return fV4.X()*std::tan(fPhiZRot/2.); } + double GetTl2() { return fV2.X()*std::tan(fPhiZRot/2.); } + + double GetH2sipm() { return fCurrentOuterHalfSipm; } + double GetBl2sipm() { return fV4sipm.X()*std::tan(fPhiZRot/2.); } + double GetTl2sipm() { return fV2sipm.X()*std::tan(fPhiZRot/2.); } + + dd4hep::RotationZYX GetRotationZYX(int numPhi); + dd4hep::Position GetTowerPos(int numPhi); + dd4hep::Position GetAssemblePos(int numPhi); + dd4hep::Position GetSipmLayerPos(int numPhi); + + dd4hep::Transform3D GetTransform3D(int numPhi); + dd4hep::Transform3D GetAssembleTransform3D(int numPhi); + dd4hep::Transform3D GetSipmTransform3D(int numPhi); + + int signedTowerNo(int unsignedTowerNo) { return fIsRHS ? unsignedTowerNo : -unsignedTowerNo-1; } + int unsignedTowerNo(int signedTowerNo) { return signedTowerNo >= 0 ? signedTowerNo : -signedTowerNo-1; } + + virtual void SetDeltaThetaByTowerNo(int , int ) {} + virtual void SetThetaOfCenterByTowerNo(int , int ) {} + void SetIsRHSByTowerNo(int signedTowerNo) { fIsRHS = ( signedTowerNo >=0 ? true : false ); } + + int GetTotTowerNum() { return fTotNum; } + void SetTotTowerNum(int totNum) { fTotNum = totNum; } + + int GetCurrentTowerNum() { return fCurrentTowerNum; } + void SetCurrentTowerNum(int numEta) { fCurrentTowerNum = numEta; } + + virtual void init() {}; + void filled() { fFilled = true; } + void finalized() { fFinalized = true; } + bool IsFinalized() { return fFinalized; } + + protected: + bool fIsRHS; + double fPhiZRot; + double fInnerX; + double fTowerH; + int fNumZRot; + double fDeltaTheta; + double fThetaOfCenter; + double fCurrentInnerR; + TVector3 fCurrentCenter; + TVector3 fV1; + TVector3 fV2; + TVector3 fV3; + TVector3 fV4; + TVector3 fV2sipm; + TVector3 fV4sipm; + double fSipmHeight; + + double fCurrentInnerHalf; + double fCurrentOuterHalf; + double fCurrentOuterHalfSipm; + + int fTotNum; + int fCurrentTowerNum; + std::vector fDeltaThetaVec; + std::vector fThetaOfCenterVec; + bool fFilled; + bool fFinalized; + }; +} +} + +#endif diff --git a/detectorSegmentations/include/detectorSegmentations/DRparamEndcap_k4geo.h b/detectorSegmentations/include/detectorSegmentations/DRparamEndcap_k4geo.h new file mode 100644 index 000000000..38572470b --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/DRparamEndcap_k4geo.h @@ -0,0 +1,24 @@ +#ifndef DETSEGMENTATION_DRPARAMENDCAP_H +#define DETSEGMENTATION_DRPARAMENDCAP_H + +#include "detectorSegmentations/DRparamBase_k4geo.h" + +#include "DD4hep/DetFactoryHelper.h" + + +namespace dd4hep { +namespace DDSegmentation { + class DRparamEndcap_k4geo : public DRparamBase_k4geo { + public: + DRparamEndcap_k4geo(); + virtual ~DRparamEndcap_k4geo(); + + virtual void SetDeltaThetaByTowerNo(int signedTowerNo, int BEtrans) override; + virtual void SetThetaOfCenterByTowerNo(int signedTowerNo, int BEtrans) override; + + virtual void init() override; + }; +} +} + +#endif diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h new file mode 100644 index 000000000..b4abef119 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h @@ -0,0 +1,164 @@ +#ifndef DETSEGMENTATION_ENDCAPTURBINE_H +#define DETSEGMENTATION_ENDCAPTURBINE_H + +// FCCSW +//#include "detectorSegmentations/GridTheta_k4geo.h" +#include "DDSegmentation/Segmentation.h" +#include "TVector3.h" + +/** FCCSWEndcapTurbine_k4geo + * + * Segmentation for turbine-style endcap calorimeter + * + */ + +namespace dd4hep { +namespace DDSegmentation { +class FCCSWEndcapTurbine_k4geo : public Segmentation { +public: + /// default constructor using an arbitrary type + FCCSWEndcapTurbine_k4geo(const std::string& aCellEncoding); + /// Default constructor used by derived classes passing an existing decoder + FCCSWEndcapTurbine_k4geo(const BitFieldCoder* decoder); + + /// destructor + virtual ~FCCSWEndcapTurbine_k4geo() = default; + + /** Determine the global position based on the cell ID. **/ + virtual Vector3D position(const CellID& aCellID) const; + /** Determine the cell ID based on the position. + * @param[in] aLocalPosition (not used). + * @param[in] aGlobalPosition position in the global coordinates. + * @param[in] aVolumeId ID of a volume. + * return Cell ID. + */ + virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, + const VolumeID& aVolumeID) const; + /** Determine the azimuthal angle based on the cell ID. + * @param[in] aCellId ID of a cell. + * return Phi. + */ + double phi(const CellID& aCellID) const; + + /** Get the coordinate offset in azimuthal angle. + * return The offset in phi. + */ + inline double offsetPhi() const { return m_offsetPhi; } + /** Get the coordinate offset in theta angle. + * return The offset in theta. + */ + inline double offsetTheta() const { return m_offsetTheta; } + /** Get the field name for azimuthal angle. + * return The field name for phi. + */ + inline const std::string& fieldNamePhi() const { return m_phiID; } + /** Set the number of bins in azimuthal angle. + * @param[in] aNumberBins Number of bins in phi. + */ + inline void setPhiBins(int bins) { m_phiBins = bins; } + /** Set the coordinate offset in azimuthal angle. + * @param[in] aOffset Offset in phi. + */ + inline void setOffsetPhi(double offset) { m_offsetPhi = offset; } + /** Set the coordinate offset in theta angle. + * @param[in] aOffset Offset in theta. + */ + inline void setOffsetTheta(double offset) { m_offsetTheta = offset; } + /** Set the field name used for phi. + * @param[in] aFieldName Field name for phi. + */ + inline void setFieldNamePhi(const std::string& fieldName) { m_phiID = fieldName; } + /** Determine the transverse distance from the beamline rho based on the cell ID. + * @param[in] aCellId ID of a cell. + * return rho. + */ + double rho(const CellID& aCellID) const; + /** Get the coordinate offset in rho. + * return The offset in rho. + */ + inline double offsetRho(int iWheel) const { return m_offsetRho[iWheel]; } + /** Get the field name for rho. + * return The field name for rho. + */ + inline const std::string& fieldNameRho() const { return m_rhoID; } + /** Set the number of bins in rho. + * @param[in] aNumberBins Number of bins in rho. + */ + inline void setRhoBins(int bins) { m_rhoBins = bins; } + /** Set the coordinate offset in rho. + * @param[in] aOffset Offset in rho. + */ + // inline void setOffsetRho(double offset) { m_offsetRho = offset; } + /** Set the field name used for transverse distance from IP rho. + * @param[in] aFieldName Field name for rho. + */ + inline void setFieldNameRho(const std::string& fieldName) { m_rhoID = fieldName; } + /** Set the field name used for the wheel ID. + * @param[in] aFieldName Field name for wheel. + */ + inline void setFieldNameWheel(const std::string& fieldName) {m_wheelID = fieldName; } + /** Determine the z coordinate based on the cell ID. + * @param[in] aCellId ID of a cell. + * return z. + */ + double z(const CellID& aCellID) const; + /** Get the coordinate offset in z. + * return The offset in z. + */ + inline double offsetZ() const { return m_offsetZ; } + /** Get the field name for z. + * return The field name for z. + */ + inline const std::string& fieldNameZ() const { return m_zID; } + /** Set the number of bins in z. + * @param[in] aNumberBins Number of bins in z. + */ + inline void setZBins(int bins) { m_zBins = bins; } + /** Set the coordinate offset in z. + * @param[in] aOffset Offset in z. + */ + inline void setOffsetZ(double offset) { m_offsetZ = offset; } + /** Set the field name used for z. + * @param[in] aFieldName Field name for z. + */ + inline void setFieldNameZ(const std::string& fieldName) { m_zID = fieldName; } + inline double rhoFromXYZ(const Vector3D& aposition) const { + TVector3 vec(aposition.X, aposition.Y, aposition.Z); + return vec.Perp(); + } + +protected: + /// turbine blade angle + double m_bladeAngle; + /// the number of bins in phi + int m_phiBins; + /// the coordinate offset in phi + double m_offsetPhi; + /// the coordinate offset in theta + double m_offsetTheta; /// the field name used for phi + std::string m_phiID; + /// the number of bins in rho + int m_rhoBins; + ////grid size in rho + std::vector m_gridSizeRho; + /// the coordinate offset in rho + std::vector m_offsetRho; + /// the field name used for rho + std::string m_rhoID; + /// the field name used for wheel + std::string m_wheelID; + /// the field name used for module + std::string m_moduleID; + /// the number of bins in z + int m_zBins; + ///grid size in z + double m_gridSizeZ; + /// the coordinate offset in z + double m_offsetZ; + /// the field name used for z + std::string m_zID; + std::string m_sideID; +}; +} +} +#endif /* DETSEGMENTATION_ENDCAPTURBINE_H */ diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h index 98fb122ca..d44f1c6c2 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h @@ -3,7 +3,6 @@ // FCCSW #include "detectorSegmentations/GridTheta_k4geo.h" -#include "DD4hep/VolumeManager.h" /** FCCSWGridModuleThetaMerged_k4geo Detector/DetSegmentation/DetSegmentation/FCCSWGridModuleThetaMerged_k4geo.h FCCSWGridModuleThetaMerged_k4geo.h * diff --git a/detectorSegmentations/include/detectorSegmentations/GridDRcaloHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridDRcaloHandle_k4geo.h new file mode 100644 index 000000000..0f777c3dd --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/GridDRcaloHandle_k4geo.h @@ -0,0 +1,90 @@ +#ifndef DD4HEP_DDCORE_GRIDDRCALO_H +#define DD4HEP_DDCORE_GRIDDRCALO_H 1 + +#include "detectorSegmentations/GridDRcalo_k4geo.h" + +#include "DD4hep/Segmentations.h" +#include "DD4hep/detail/SegmentationsInterna.h" + +namespace dd4hep { + class Segmentation; + template class SegmentationWrapper; + + typedef Handle> GridDRcaloHandle; + + class GridDRcalo_k4geo : public GridDRcaloHandle { + typedef GridDRcaloHandle::Object Object; + + public: + /// Default constructor + GridDRcalo_k4geo() = default; + /// Copy constructor + GridDRcalo_k4geo(const GridDRcalo_k4geo& e) = default; + /// Copy Constructor from segmentation base object + GridDRcalo_k4geo(const Segmentation& e) : Handle(e) {} + /// Copy constructor from handle + GridDRcalo_k4geo(const Handle& e) : Handle(e) {} + /// Copy constructor from other polymorph/equivalent handle + template + GridDRcalo_k4geo(const Handle& e) : Handle(e) {} + /// Assignment operator + GridDRcalo_k4geo& operator=(const GridDRcalo_k4geo& seg) = default; + /// Equality operator + bool operator==(const GridDRcalo_k4geo& seg) const { return m_element == seg.m_element; } + /// determine the position based on the cell ID + inline Position position(const CellID& id) const { return Position(access()->implementation->position(id)); } + inline Position localPosition(const CellID& id) const { return Position(access()->implementation->localPosition(id)); } + inline Position localPosition(int numx, int numy, int x_, int y_) const { return Position(access()->implementation->localPosition(numx,numy,x_,y_)); } + + /// determine the cell ID based on the position + inline dd4hep::CellID cellID(const Position& local, const Position& global, const VolumeID& volID) const { + return access()->implementation->cellID(local, global, volID); + } + + inline VolumeID setVolumeID(int numEta, int numPhi) const { return access()->implementation->setVolumeID(numEta,numPhi); } + inline CellID setCellID(int numEta, int numPhi, int x, int y) const { return access()->implementation->setCellID(numEta, numPhi, x, y); } + + inline void setGridSize(double grid) { access()->implementation->setGridSize(grid); } + + // Get the identifier number of a mother tower in eta or phi direction + inline int numEta(const CellID& aCellID) const { return access()->implementation->numEta(aCellID); } + inline int numPhi(const CellID& aCellID) const { return access()->implementation->numPhi(aCellID); } + + // Get the total number of SiPMs of the mother tower in x or y direction (local coordinate) + inline int numX(const CellID& aCellID) const { return access()->implementation->numX(aCellID); } + inline int numY(const CellID& aCellID) const { return access()->implementation->numY(aCellID); } + + // Get the identifier number of a SiPM in x or y direction (local coordinate) + inline int x(const CellID& aCellID) const { return access()->implementation->x(aCellID); } // approx eta direction + inline int y(const CellID& aCellID) const { return access()->implementation->y(aCellID); } // approx phi direction + + inline bool IsCerenkov(const CellID& aCellID) const { return access()->implementation->IsCerenkov(aCellID); } + inline bool IsCerenkov(int col, int row) const { return access()->implementation->IsCerenkov(col, row); } + + inline bool IsTower(const CellID& aCellID) const { return access()->implementation->IsTower(aCellID); } + inline bool IsSiPM(const CellID& aCellID) const { return access()->implementation->IsSiPM(aCellID); } + + inline int getFirst32bits(const CellID& aCellID) const { return access()->implementation->getFirst32bits(aCellID); } + inline int getLast32bits(const CellID& aCellID) const { return access()->implementation->getLast32bits(aCellID); } + + inline CellID convertFirst32to64(const int aId32) const { return access()->implementation->convertFirst32to64(aId32); } + inline CellID convertLast32to64(const int aId32) const { return access()->implementation->convertLast32to64(aId32); } + + // Methods for 32bit to 64bit en/decoder + inline int numEta(const int& aId32) const { return access()->implementation->numEta(aId32); } + inline int numPhi(const int& aId32) const { return access()->implementation->numPhi(aId32); } + + inline int numX(const int& aId32) const { return access()->implementation->numX(aId32); } + inline int numY(const int& aId32) const { return access()->implementation->numY(aId32); } + + inline int x(const int& aId32) const { return access()->implementation->x(aId32); } + inline int y(const int& aId32) const { return access()->implementation->y(aId32); } + + inline bool IsCerenkov(const int& aId32) const { return access()->implementation->IsCerenkov(aId32); } + + inline bool IsTower(const int& aId32) const { return access()->implementation->IsTower(aId32); } + inline bool IsSiPM(const int& aId32) const { return access()->implementation->IsSiPM(aId32); } + }; +} + +#endif diff --git a/detectorSegmentations/include/detectorSegmentations/GridDRcalo_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridDRcalo_k4geo.h new file mode 100644 index 000000000..52eb530e5 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/GridDRcalo_k4geo.h @@ -0,0 +1,112 @@ +#ifndef DETSEGMENTATION_GRIDDRCALO_H +#define DETSEGMENTATION_GRIDDRCALO_H + +#include "detectorSegmentations/DRparamBarrel_k4geo.h" +#include "detectorSegmentations/DRparamEndcap_k4geo.h" + +#include "DDSegmentation/Segmentation.h" + +namespace dd4hep { +namespace DDSegmentation { +class GridDRcalo_k4geo : public Segmentation { +public: + /// default constructor using an arbitrary type + GridDRcalo_k4geo(const std::string& aCellEncoding); + /// Default constructor used by derived classes passing an existing decoder + GridDRcalo_k4geo(const BitFieldCoder* decoder); + /// destructor + virtual ~GridDRcalo_k4geo() override; + + // Determine the global(local) position based on the cell ID. + virtual Vector3D position(const CellID& aCellID) const override; + Vector3D localPosition(const CellID& aCellID) const; + Vector3D localPosition(int numx, int numy, int x_, int y_) const; + + virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, + const VolumeID& aVolumeID) const override; + + VolumeID setVolumeID(int numEta, int numPhi) const; + CellID setCellID(int numEta, int numPhi, int x, int y) const; + + void setGridSize(double grid) { fGridSize = grid; } + void setSipmSize(double sipm) { fSipmSize = sipm; } + + // Get the identifier number of a mother tower in eta or phi direction + int numEta(const CellID& aCellID) const; + int numPhi(const CellID& aCellID) const; + + // Get the total number of SiPMs of the mother tower in x or y direction (local coordinate) + int numX(const CellID& aCellID) const; + int numY(const CellID& aCellID) const; + + // Get the identifier number of a SiPM in x or y direction (local coordinate) + int x(const CellID& aCellID) const; // approx eta direction + int y(const CellID& aCellID) const; // approx phi direction + + bool IsCerenkov(const CellID& aCellID) const; + bool IsCerenkov(int col, int row) const; + + bool IsTower(const CellID& aCellID) const; + bool IsSiPM(const CellID& aCellID) const; + + int getFirst32bits(const CellID& aCellID) const { return (int)aCellID; } + int getLast32bits(const CellID& aCellID) const; + CellID convertFirst32to64(const int aId32) const { return (CellID)aId32; } + CellID convertLast32to64(const int aId32) const; + + // Methods for 32bit to 64bit en/decoder + int numEta(const int& aId32) const { return numEta( convertFirst32to64(aId32) ); } + int numPhi(const int& aId32) const { return numPhi( convertFirst32to64(aId32) ); } + + int numX(const int& aId32) const { return numX( convertFirst32to64(aId32) ); } + int numY(const int& aId32) const { return numY( convertFirst32to64(aId32) ); } + + int x(const int& aId32) const { return x( convertLast32to64(aId32) ); } + int y(const int& aId32) const { return y( convertLast32to64(aId32) ); } + + bool IsCerenkov(const int& aId32) const { return IsCerenkov( convertLast32to64(aId32) ); } + + bool IsTower(const int& aId32) const { return IsTower( convertLast32to64(aId32) ); } + bool IsSiPM(const int& aId32) const { return IsSiPM( convertLast32to64(aId32) ); } + + inline const std::string& fieldNameAssembly() const { return fAssemblyId; } + inline const std::string& fieldNameNumEta() const { return fNumEtaId; } + inline const std::string& fieldNameNumPhi() const { return fNumPhiId; } + inline const std::string& fieldNameX() const { return fXId; } + inline const std::string& fieldNameY() const { return fYId; } + inline const std::string& fieldNameIsCerenkov() const { return fIsCerenkovId; } + inline const std::string& fieldNameModule() const { return fModule; } + + inline void setFieldNameAssembly(const std::string& fieldName) { fAssemblyId = fieldName; } + inline void setFieldNameNumEta(const std::string& fieldName) { fNumEtaId = fieldName; } + inline void setFieldNameNumPhi(const std::string& fieldName) { fNumPhiId = fieldName; } + inline void setFieldNameX(const std::string& fieldName) { fXId = fieldName; } + inline void setFieldNameY(const std::string& fieldName) { fYId = fieldName; } + inline void setFieldNameIsCerenkov(const std::string& fieldName) { fIsCerenkovId = fieldName; } + inline void setFieldNameModule(const std::string& fieldName) { fModule = fieldName; } + + DRparamBarrel_k4geo* paramBarrel() { return fParamBarrel; } + DRparamEndcap_k4geo* paramEndcap() { return fParamEndcap; } + + DRparamBase_k4geo* setParamBase(int noEta) const; + +protected: + std::string fAssemblyId; + std::string fNumEtaId; + std::string fNumPhiId; + std::string fXId; + std::string fYId; + std::string fIsCerenkovId; + std::string fModule; + + double fGridSize; + double fSipmSize; + +private: + DRparamBarrel_k4geo* fParamBarrel; + DRparamEndcap_k4geo* fParamEndcap; +}; +} +} + +#endif diff --git a/detectorSegmentations/include/detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h index dcb389b1c..be1b43f6d 100644 --- a/detectorSegmentations/include/detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h @@ -5,7 +5,6 @@ #include "TVector3.h" #include -#include /** GridSimplifiedDriftChamber_k4geo Detector/detectorSegmentations/detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h GridSimplifiedDriftChamber_k4geo.h * diff --git a/detectorSegmentations/src/DRparamBarrel_k4geo.cpp b/detectorSegmentations/src/DRparamBarrel_k4geo.cpp new file mode 100644 index 000000000..3daf8e681 --- /dev/null +++ b/detectorSegmentations/src/DRparamBarrel_k4geo.cpp @@ -0,0 +1,97 @@ +#include "detectorSegmentations/DRparamBarrel_k4geo.h" + + +#include + +namespace dd4hep { +namespace DDSegmentation { + +DRparamBarrel_k4geo::DRparamBarrel_k4geo() { + fIsRHS = 0; + fPhiZRot = 0.; + fInnerX = 0.; + fTowerH = 0.; + fNumZRot = 0; + fDeltaTheta = 0.; + fThetaOfCenter = 0.; + fCurrentInnerR = 0.; + fPhiZRot = 0; + fCurrentCenter = TVector3(); + fV1 = TVector3(); + fV2 = TVector3(); + fV3 = TVector3(); + fV4 = TVector3(); + fSipmHeight = 0.; + fCurrentInnerHalf = 0.; + fCurrentOuterHalf = 0.; + fFilled = false; + fFinalized = false; +} + +DRparamBarrel_k4geo::~DRparamBarrel_k4geo() {} + +void DRparamBarrel_k4geo::init() { + fCurrentInnerR = fInnerX/std::cos(fThetaOfCenter); + double trnsLength = fTowerH/2.+fCurrentInnerR; + fCurrentCenter = TVector3(std::cos(fThetaOfCenter)*trnsLength,0.,std::sin(fThetaOfCenter)*trnsLength); + + fCurrentInnerHalf = fCurrentInnerR*std::tan(fDeltaTheta/2.); + fCurrentOuterHalf = (fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.); + fCurrentOuterHalfSipm = (fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.); + + fV1 = TVector3( + std::cos(fThetaOfCenter)*fCurrentInnerR+std::sin(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*fCurrentInnerR-std::cos(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.) + ); + + fV2 = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)+std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)-std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.) + ); + + fV3 = TVector3( + std::cos(fThetaOfCenter)*fCurrentInnerR-std::sin(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*fCurrentInnerR+std::cos(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.) + ); + + fV4 = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)-std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)+std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.) + ); + + fV2sipm = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)+std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)-std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.) + ); + + fV4sipm = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)-std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)+std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.) + ); + + if (!fFilled) { + fDeltaThetaVec.push_back(fDeltaTheta); + fThetaOfCenterVec.push_back(fThetaOfCenter); + } +} + +void DRparamBarrel_k4geo::SetDeltaThetaByTowerNo(int signedTowerNo, int) { + if (!fFilled) throw std::runtime_error("Attempt to set by tower num while barrel parameter is not filled!"); + + fDeltaTheta = fDeltaThetaVec.at( unsignedTowerNo(signedTowerNo) ); +} + +void DRparamBarrel_k4geo::SetThetaOfCenterByTowerNo(int signedTowerNo, int) { + if (!fFilled) throw std::runtime_error("Attempt to set by tower num while barrel parameter is not filled!"); + + fThetaOfCenter = fThetaOfCenterVec.at( unsignedTowerNo(signedTowerNo) ); +} + +} +} diff --git a/detectorSegmentations/src/DRparamBase_k4geo.cpp b/detectorSegmentations/src/DRparamBase_k4geo.cpp new file mode 100644 index 000000000..964362947 --- /dev/null +++ b/detectorSegmentations/src/DRparamBase_k4geo.cpp @@ -0,0 +1,99 @@ +#include "detectorSegmentations/DRparamBase_k4geo.h" + +#include "Math/GenVector/RotationZYX.h" + + +namespace dd4hep { +namespace DDSegmentation { + +DRparamBase_k4geo::DRparamBase_k4geo() { + fIsRHS = 0; + fPhiZRot = 0.; + fInnerX = 0.; + fTowerH = 0.; + fNumZRot = 0; + fDeltaTheta = 0.; + fThetaOfCenter = 0.; + fCurrentInnerR = 0.; + fPhiZRot = 0; + fCurrentCenter = TVector3(); + fV1 = TVector3(); + fV2 = TVector3(); + fV3 = TVector3(); + fV4 = TVector3(); + fSipmHeight = 0.; + fCurrentInnerHalf = 0.; + fCurrentOuterHalf = 0.; + fCurrentTowerNum = 0; + fFilled = false; + fFinalized = false; +} + +DRparamBase_k4geo::~DRparamBase_k4geo() {} + +dd4hep::RotationZYX DRparamBase_k4geo::GetRotationZYX(int numPhi) { + double numPhi_ = (double)numPhi; + double xRot = fIsRHS ? -fThetaOfCenter : fThetaOfCenter; + double zRot = fIsRHS ? -M_PI/2. : M_PI/2.; + dd4hep::RotationZYX rot = dd4hep::RotationZYX(zRot, M_PI/2.+xRot, 0.); + ROOT::Math::RotationZ rotZ = ROOT::Math::RotationZ(numPhi_*fPhiZRot); + rot = rotZ*rot; + + return rot; +} + +dd4hep::Position DRparamBase_k4geo::GetTowerPos(int numPhi) { + double numPhi_ = (double)numPhi; + double x = std::cos(numPhi_*fPhiZRot)*fCurrentCenter.X(); + double y = std::sin(numPhi_*fPhiZRot)*fCurrentCenter.X(); + double z = fIsRHS ? fCurrentCenter.Z() : -fCurrentCenter.Z(); + dd4hep::Position pos = dd4hep::Position(x,y,z); + + return pos; +} + +dd4hep::Position DRparamBase_k4geo::GetAssemblePos(int numPhi) { + double numPhi_ = (double)numPhi; + double x = std::cos(numPhi_*fPhiZRot)*fCurrentCenter.X()*(fCurrentCenter.Mag()+fSipmHeight/2.)/fCurrentCenter.Mag(); + double y = std::sin(numPhi_*fPhiZRot)*fCurrentCenter.X()*(fCurrentCenter.Mag()+fSipmHeight/2.)/fCurrentCenter.Mag(); + double z_abs = fCurrentCenter.Z()*(fCurrentCenter.Mag()+fSipmHeight/2.)/fCurrentCenter.Mag(); + double z = fIsRHS ? z_abs : -z_abs; + dd4hep::Position pos = dd4hep::Position(x,y,z); + + return pos; +} + +dd4hep::Position DRparamBase_k4geo::GetSipmLayerPos(int numPhi) { + double numPhi_ = (double)numPhi; + double x = std::cos(numPhi_*fPhiZRot)*fCurrentCenter.X()*(fCurrentCenter.Mag()+fTowerH/2.+fSipmHeight/2.)/fCurrentCenter.Mag(); + double y = std::sin(numPhi_*fPhiZRot)*fCurrentCenter.X()*(fCurrentCenter.Mag()+fTowerH/2.+fSipmHeight/2.)/fCurrentCenter.Mag(); + double z_abs = fCurrentCenter.Z()*(fCurrentCenter.Mag()+fTowerH/2.+fSipmHeight/2.)/fCurrentCenter.Mag(); + double z = fIsRHS ? z_abs : -z_abs; + dd4hep::Position pos = dd4hep::Position(x,y,z); + + return pos; +} + +dd4hep::Transform3D DRparamBase_k4geo::GetTransform3D(int numPhi) { + auto rot = GetRotationZYX(numPhi); + auto pos = GetTowerPos(numPhi); + + return dd4hep::Transform3D(rot,pos); +} + +dd4hep::Transform3D DRparamBase_k4geo::GetAssembleTransform3D(int numPhi) { + auto rot = GetRotationZYX(numPhi); + auto pos = GetAssemblePos(numPhi); + + return dd4hep::Transform3D(rot,pos); +} + +dd4hep::Transform3D DRparamBase_k4geo::GetSipmTransform3D(int numPhi) { + auto rot = GetRotationZYX(numPhi); + auto pos = GetSipmLayerPos(numPhi); + + return dd4hep::Transform3D(rot,pos); +} + +} +} diff --git a/detectorSegmentations/src/DRparamEndcap_k4geo.cpp b/detectorSegmentations/src/DRparamEndcap_k4geo.cpp new file mode 100644 index 000000000..09f9b78c5 --- /dev/null +++ b/detectorSegmentations/src/DRparamEndcap_k4geo.cpp @@ -0,0 +1,97 @@ +#include "detectorSegmentations/DRparamEndcap_k4geo.h" + + +#include + +namespace dd4hep { +namespace DDSegmentation { + +DRparamEndcap_k4geo::DRparamEndcap_k4geo() { + fIsRHS = 0; + fPhiZRot = 0.; + fInnerX = 0.; + fTowerH = 0.; + fNumZRot = 0; + fDeltaTheta = 0.; + fThetaOfCenter = 0.; + fCurrentInnerR = 0.; + fPhiZRot = 0; + fCurrentCenter = TVector3(); + fV1 = TVector3(); + fV2 = TVector3(); + fV3 = TVector3(); + fV4 = TVector3(); + fSipmHeight = 0.; + fCurrentInnerHalf = 0.; + fCurrentOuterHalf = 0.; + fFilled = false; + fFinalized = false; +} + +DRparamEndcap_k4geo::~DRparamEndcap_k4geo() {} + +void DRparamEndcap_k4geo::init() { + fCurrentInnerR = fInnerX/std::sin(fThetaOfCenter); + double trnsLength = fTowerH/2.+fCurrentInnerR; + fCurrentCenter = TVector3(std::cos(fThetaOfCenter)*trnsLength,0.,std::sin(fThetaOfCenter)*trnsLength); + + fCurrentInnerHalf = fCurrentInnerR*std::tan(fDeltaTheta/2.); + fCurrentOuterHalf = (fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.); + fCurrentOuterHalfSipm = (fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.); + + fV1 = TVector3( + std::cos(fThetaOfCenter)*fCurrentInnerR+std::sin(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*fCurrentInnerR-std::cos(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.) + ); + + fV2 = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)+std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)-std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.) + ); + + fV3 = TVector3( + std::cos(fThetaOfCenter)*fCurrentInnerR-std::sin(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*fCurrentInnerR+std::cos(fThetaOfCenter)*fCurrentInnerR*std::tan(fDeltaTheta/2.) + ); + + fV4 = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)-std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH)+std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH)*std::tan(fDeltaTheta/2.) + ); + + fV2sipm = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)+std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)-std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.) + ); + + fV4sipm = TVector3( + std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)-std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.), + 0., + std::sin(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)+std::cos(fThetaOfCenter)*(fCurrentInnerR+fTowerH+fSipmHeight)*std::tan(fDeltaTheta/2.) + ); + + if (!fFilled) { + fDeltaThetaVec.push_back(fDeltaTheta); + fThetaOfCenterVec.push_back(fThetaOfCenter); + } +} + +void DRparamEndcap_k4geo::SetDeltaThetaByTowerNo(int signedTowerNo, int BEtrans) { + if (!fFilled) throw std::runtime_error("Attempt to set by tower num while endcap parameter is not filled!"); + + fDeltaTheta = fDeltaThetaVec.at( unsignedTowerNo(signedTowerNo) - unsignedTowerNo(BEtrans) ); +} + +void DRparamEndcap_k4geo::SetThetaOfCenterByTowerNo(int signedTowerNo, int BEtrans) { + if (!fFilled) throw std::runtime_error("Attempt to set by tower num while endcap parameter is not filled!"); + + fThetaOfCenter = fThetaOfCenterVec.at( unsignedTowerNo(signedTowerNo) - unsignedTowerNo(BEtrans) ); +} + +} +} diff --git a/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp b/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp new file mode 100644 index 000000000..c25d00625 --- /dev/null +++ b/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp @@ -0,0 +1,123 @@ +#include "detectorSegmentations/FCCSWEndcapTurbine_k4geo.h" +#include "DD4hep/Detector.h" + +namespace dd4hep { +namespace DDSegmentation { + +/// default constructor using an encoding string + FCCSWEndcapTurbine_k4geo::FCCSWEndcapTurbine_k4geo(const std::string& cellEncoding) : Segmentation(cellEncoding) { + // define type and description + _type = "FCCSWEndcapTurbine_k4geo"; + _description = "Turbine-specific segmentation in the global coordinates"; + + // register all necessary parameters + registerParameter("grid_size_rho", "Grid size in rho", m_gridSizeRho, std::vector()); + registerParameter("offset_rho", "Offset in rho", m_offsetRho, std::vector()); + registerIdentifier("identifier_rho", "Cell ID identifier for rho", m_rhoID, "rho"); // might want to have separate concepts for rho and layer... + registerParameter("grid_size_z", "Grid size in z", m_gridSizeZ, 0.); + registerParameter("offset_z", "Offset in z", m_offsetZ, 0.); + registerParameter("offset_theta", "Angular offset in theta", m_offsetTheta, 0., SegmentationParameter::AngleUnit, true); + registerIdentifier("identifier_z", "Cell ID identifier for z", m_zID, "z"); + registerIdentifier("identifier_side", "Cell ID identifier for side", m_sideID, "side"); + + dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); + try { + m_bladeAngle = dd4hepgeo->constant("BladeAngle"); + std::cout << "Blade angle is " << m_bladeAngle << std::endl; + } + catch(...) { + std::cout << "Blade angle not found in detector metadata, exiting..." << std::endl; + exit(1); + } + +} + + FCCSWEndcapTurbine_k4geo::FCCSWEndcapTurbine_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { + // define type and description + _type = "FCCSWEndcapTurbine_k4geo"; + _description = "Turbine-specific segmentation in the global coordinates"; + + // register all necessary parameters + registerParameter("grid_size_rho", "Grid size in rho", m_gridSizeRho, std::vector()); + registerParameter("offset_rho", "Offset in rho", m_offsetRho, std::vector()); + registerIdentifier("identifier_rho", "Cell ID identifier for rho", m_rhoID, "rho"); + registerParameter("grid_size_z", "Grid size in z", m_gridSizeZ, 0.); + registerParameter("offset_z", "Offset in z", m_offsetZ, 0.); + registerParameter("offset_theta", "Angular offset in theta", m_offsetTheta, 0., SegmentationParameter::AngleUnit, true); + registerIdentifier("identifier_z", "Cell ID identifier for z", m_zID, "z"); + registerIdentifier("identifier_side", "Cell ID identifier for side", m_sideID, "side"); + registerIdentifier("identifier_wheel", "Cell ID identifier for wheel", m_wheelID, "wheel"); + registerIdentifier("identifier_module", "Cell ID identifier for module", m_moduleID, "module"); + dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); + try { + m_bladeAngle = dd4hepgeo->constant("BladeAngle"); + std::cout << "Blade angle is " << m_bladeAngle << std::endl; + } + catch(...) { + std::cout << "Blade angle not found in detector metadata, exiting..." << std::endl; + exit(1); + } + + +} + +/// determine the local position based on the cell ID +Vector3D FCCSWEndcapTurbine_k4geo::position(const CellID& cID) const { + + double rhoVal = rho(cID); + double zVal = z(cID); + double phiVal = phi(cID); + Vector3D pos = PositionRhoZPhi(rhoVal, zVal, phiVal); + // account for the fact that the -z endcap is mirrored wrt to the +z one + if (pos.Z < 0.) pos.Y = -pos.Y; + return pos; +} + +/// determine the cell ID based on the position +CellID FCCSWEndcapTurbine_k4geo::cellID(const Vector3D& /* localPosition */, const Vector3D& globalPosition, + const VolumeID& vID) const { + CellID cID = vID; + double lRho = rhoFromXYZ(globalPosition); + CellID iWheel = _decoder->get(cID, m_wheelID); + _decoder->set(cID, m_rhoID, positionToBin(lRho, m_gridSizeRho[iWheel], m_offsetRho[iWheel])); + double lZ = TMath::Abs(globalPosition.Z); + _decoder->set(cID, m_zID, positionToBin(lZ, m_gridSizeZ, m_offsetZ)); + + return cID; +} + +/// determine the azimuthal angle phi based on the cell ID +double FCCSWEndcapTurbine_k4geo::phi(const CellID& cID) const { + + CellID iModule = _decoder->get(cID, m_moduleID); + double phiCent = twopi*(iModule+0.5)/78336; + + double rhoLoc = rho(cID); + double zLoc = TMath::Abs(z(cID))-m_offsetZ - 45; // hard-code midpoint in z for now + + double zCotBladeAngle = zLoc/TMath::Tan(m_bladeAngle); + double x = zCotBladeAngle; + double y = TMath::Sqrt(rhoLoc*rhoLoc - x*x); + // rotate about z axis by phiCent + double xprime = x*TMath::Cos(phiCent) +y*TMath::Sin(phiCent); + double yprime = y*TMath::Cos(phiCent) -x*TMath::Sin(phiCent); + + return TMath::ATan2(xprime,yprime); +} +/// determine the transverse distance from the beamline r based on the cell ID +double FCCSWEndcapTurbine_k4geo::rho(const CellID& cID) const { + CellID rhoValue = _decoder->get(cID, m_rhoID); + CellID iWheel = _decoder->get(cID, m_wheelID); + + return binToPosition(rhoValue,m_gridSizeRho[iWheel], m_offsetRho[iWheel]); +} +/// determine local x in plane of blade based on the cell ID +double FCCSWEndcapTurbine_k4geo::z(const CellID& cID) const { + CellID zValue = _decoder->get(cID, m_zID); + CellID sideValue = _decoder->get(cID, m_sideID); + return ((long long int)sideValue)*binToPosition(zValue,m_gridSizeZ,m_offsetZ); +} + +} +} + diff --git a/detectorSegmentations/src/GridDRcaloHandle_k4geo.cpp b/detectorSegmentations/src/GridDRcaloHandle_k4geo.cpp new file mode 100644 index 000000000..88bc40586 --- /dev/null +++ b/detectorSegmentations/src/GridDRcaloHandle_k4geo.cpp @@ -0,0 +1,4 @@ +#include "detectorSegmentations/GridDRcaloHandle_k4geo.h" +#include "DD4hep/detail/Handle.inl" + +DD4HEP_INSTANTIATE_HANDLE_UNNAMED(dd4hep::DDSegmentation::GridDRcalo_k4geo); diff --git a/detectorSegmentations/src/GridDRcalo_k4geo.cpp b/detectorSegmentations/src/GridDRcalo_k4geo.cpp new file mode 100644 index 000000000..ef4cc90e6 --- /dev/null +++ b/detectorSegmentations/src/GridDRcalo_k4geo.cpp @@ -0,0 +1,239 @@ +#include "detectorSegmentations/GridDRcalo_k4geo.h" + +#include +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + +/// default constructor using an encoding string +GridDRcalo_k4geo::GridDRcalo_k4geo(const std::string& cellEncoding) : Segmentation(cellEncoding) { + // define type and description + _type = "GridDRcalo_k4geo"; + _description = "DRcalo segmentation based on the tower / (Cerenkov or Scintillation) fiber / SiPM hierarchy"; + + // register all necessary parameters + registerIdentifier("identifier_assembly", "Cell ID identifier for assembly", fAssemblyId, "assembly"); + registerIdentifier("identifier_eta", "Cell ID identifier for numEta", fNumEtaId, "eta"); + registerIdentifier("identifier_phi", "Cell ID identifier for numPhi", fNumPhiId, "phi"); + registerIdentifier("identifier_x", "Cell ID identifier for x", fXId, "x"); + registerIdentifier("identifier_y", "Cell ID identifier for y", fYId, "y"); + registerIdentifier("identifier_IsCerenkov", "Cell ID identifier for IsCerenkov", fIsCerenkovId, "c"); + registerIdentifier("identifier_module", "Cell ID identifier for module", fModule, "module"); + + fParamBarrel = new DRparamBarrel_k4geo(); + fParamEndcap = new DRparamEndcap_k4geo(); +} + +GridDRcalo_k4geo::GridDRcalo_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { + // define type and description + _type = "GridDRcalo_k4geo"; + _description = "DRcalo segmentation based on the tower / (Cerenkov or Scintillation) fiber / SiPM hierarchy"; + + // register all necessary parameters + registerIdentifier("identifier_assembly", "Cell ID identifier for assembly", fAssemblyId, "assembly"); + registerIdentifier("identifier_eta", "Cell ID identifier for numEta", fNumEtaId, "eta"); + registerIdentifier("identifier_phi", "Cell ID identifier for numPhi", fNumPhiId, "phi"); + registerIdentifier("identifier_x", "Cell ID identifier for x", fXId, "x"); + registerIdentifier("identifier_y", "Cell ID identifier for y", fYId, "y"); + registerIdentifier("identifier_IsCerenkov", "Cell ID identifier for IsCerenkov", fIsCerenkovId, "c"); + registerIdentifier("identifier_module", "Cell ID identifier for module", fModule, "module"); + + fParamBarrel = new DRparamBarrel_k4geo(); + fParamEndcap = new DRparamEndcap_k4geo(); +} + +GridDRcalo_k4geo::~GridDRcalo_k4geo() { + delete fParamBarrel; + delete fParamEndcap; +} + +Vector3D GridDRcalo_k4geo::position(const CellID& cID) const { + int noEta = numEta(cID); + int noPhi = numPhi(cID); + + DRparamBase_k4geo* paramBase = setParamBase(noEta); + + auto transformA = paramBase->GetSipmTransform3D(noPhi); + dd4hep::Position localPos = dd4hep::Position(0.,0.,0.); + if ( IsSiPM(cID) ) localPos = dd4hep::Position( localPosition(cID) ); + + dd4hep::RotationZYX rot = dd4hep::RotationZYX(M_PI, 0., 0.); // AdHoc rotation, potentially bug + dd4hep::Transform3D transformB = dd4hep::Transform3D(rot,localPos); + auto total = transformA*transformB; + + dd4hep::Position translation = dd4hep::Position(0.,0.,0.); + total.GetTranslation(translation); + + return Vector3D(translation.x(),translation.y(),translation.z()); +} + +Vector3D GridDRcalo_k4geo::localPosition(const CellID& cID) const { + int numx = numX(cID); + int numy = numY(cID); + int x_ = x(cID); + int y_ = y(cID); + + return localPosition(numx,numy,x_,y_); +} + +Vector3D GridDRcalo_k4geo::localPosition(int numx, int numy, int x_, int y_) const { + // Here x & y range from 0 to numx & numy (used for the geometry construction) + float ptX = -fGridSize*static_cast(numx/2) + static_cast(x_)*fGridSize + ( numx%2==0 ? fGridSize/2. : 0. ); + float ptY = -fGridSize*static_cast(numy/2) + static_cast(y_)*fGridSize + ( numy%2==0 ? fGridSize/2. : 0. ); + + return Vector3D(ptX,ptY,0.); +} + +/// determine the cell ID based on the position +CellID GridDRcalo_k4geo::cellID(const Vector3D& localPosition, const Vector3D& globalPosition, const VolumeID& vID) const { + int numx = numX(vID); + int numy = numY(vID); + + auto localX = localPosition.x(); + auto localY = localPosition.y(); + + int x = std::floor( ( localX + ( numx%2==0 ? 0. : fGridSize/2. ) ) / fGridSize ) + numx/2; + int y = std::floor( ( localY + ( numy%2==0 ? 0. : fGridSize/2. ) ) / fGridSize ) + numy/2; + + int signedNumEta = (globalPosition.z() >= 0) ? numEta(vID) : (-numEta(vID) - 1); + return setCellID( signedNumEta, numPhi(vID), x, y ); +} + +VolumeID GridDRcalo_k4geo::setVolumeID(int numEta, int numPhi) const { + VolumeID numEtaId = static_cast(numEta); + VolumeID numPhiId = static_cast(numPhi); + VolumeID vID = 0; + _decoder->set(vID, fNumEtaId, numEtaId); + _decoder->set(vID, fNumPhiId, numPhiId); + + VolumeID module = 0; // Tower, SiPM layer attached to the tower, etc. + _decoder->set(vID, fModule, module); + + return vID; +} + +CellID GridDRcalo_k4geo::setCellID(int numEta, int numPhi, int x, int y) const { + VolumeID numEtaId = static_cast(numEta); + VolumeID numPhiId = static_cast(numPhi); + VolumeID xId = static_cast(x); + VolumeID yId = static_cast(y); + VolumeID vID = 0; + _decoder->set(vID, fNumEtaId, numEtaId); + _decoder->set(vID, fNumPhiId, numPhiId); + _decoder->set(vID, fXId, xId); + _decoder->set(vID, fYId, yId); + + VolumeID module = 1; // Fiber, SiPM, etc. + _decoder->set(vID, fModule, module); + + VolumeID isCeren = IsCerenkov(x,y) ? 1 : 0; + _decoder->set(vID, fIsCerenkovId, isCeren); + + VolumeID assemblyId = (numEta >= 0) ? 0 : 1; + _decoder->set(vID, fAssemblyId, assemblyId); + + return vID; +} + +// Get the identifier number of a mother tower in eta or phi direction +int GridDRcalo_k4geo::numEta(const CellID& aCellID) const { + VolumeID numEta = static_cast(_decoder->get(aCellID, fNumEtaId)); + return static_cast(numEta); +} + +int GridDRcalo_k4geo::numPhi(const CellID& aCellID) const { + VolumeID numPhi = static_cast(_decoder->get(aCellID, fNumPhiId)); + return static_cast(numPhi); +} + +// Get the total number of SiPMs of the mother tower in x or y direction (local coordinate) +int GridDRcalo_k4geo::numX(const CellID& aCellID) const { + int noEta = numEta(aCellID); + + DRparamBase_k4geo* paramBase = setParamBase(noEta); + + int noX = static_cast( std::floor( ( paramBase->GetTl2()*2. - fSipmSize/2. )/fGridSize ) ) + 1; // in phi direction + + return noX; +} + +int GridDRcalo_k4geo::numY(const CellID& aCellID) const { + int noEta = numEta(aCellID); + + DRparamBase_k4geo* paramBase = setParamBase(noEta); + + int noY = static_cast( std::floor( ( paramBase->GetH2()*2. - fSipmSize/2. )/fGridSize ) ) + 1; // in eta direction + + return noY; +} + +// Get the identifier number of a SiPM in x or y direction (local coordinate) +int GridDRcalo_k4geo::x(const CellID& aCellID) const { // approx eta direction + VolumeID x = static_cast(_decoder->get(aCellID, fXId)); + return static_cast(x); +} +int GridDRcalo_k4geo::y(const CellID& aCellID) const { // approx phi direction + VolumeID y = static_cast(_decoder->get(aCellID, fYId)); + return static_cast(y); +} + +bool GridDRcalo_k4geo::IsCerenkov(const CellID& aCellID) const { + VolumeID isCeren = static_cast(_decoder->get(aCellID, fIsCerenkovId)); + return static_cast(isCeren); +} +// Identify if the fiber is Cerenkov or scintillation by its column and row number +bool GridDRcalo_k4geo::IsCerenkov(int col, int row) const { + bool isCeren = false; + if ( col%2 == 1 ) { isCeren = !isCeren; } + if ( row%2 == 1 ) { isCeren = !isCeren; } + return isCeren; +} + +bool GridDRcalo_k4geo::IsTower(const CellID& aCellID) const { + VolumeID module = static_cast(_decoder->get(aCellID, fModule)); + return module==0; +} + +bool GridDRcalo_k4geo::IsSiPM(const CellID& aCellID) const { + VolumeID module = static_cast(_decoder->get(aCellID, fModule)); + return module==1; +} + +int GridDRcalo_k4geo::getLast32bits(const CellID& aCellID) const { + CellID aId64 = aCellID >> sizeof(int)*CHAR_BIT; + int aId32 = (int)aId64; + + return aId32; +} + +CellID GridDRcalo_k4geo::convertLast32to64(const int aId32) const { + CellID aId64 = (CellID)aId32; + aId64 <<= sizeof(int)*CHAR_BIT; + + return aId64; +} + +DRparamBase_k4geo* GridDRcalo_k4geo::setParamBase(int noEta) const { + DRparamBase_k4geo* paramBase = nullptr; + + if ( fParamEndcap->unsignedTowerNo(noEta) >= fParamBarrel->GetTotTowerNum() ) paramBase = static_cast(fParamEndcap); + else paramBase = static_cast(fParamBarrel); + + if ( paramBase->GetCurrentTowerNum()==noEta ) return paramBase; + + // This should not be called while building detector geometry + if (!paramBase->IsFinalized()) throw std::runtime_error("GridDRcalo_k4geo::position should not be called while building detector geometry!"); + + paramBase->SetDeltaThetaByTowerNo(noEta, fParamBarrel->GetTotTowerNum()); + paramBase->SetThetaOfCenterByTowerNo(noEta, fParamBarrel->GetTotTowerNum()); + paramBase->SetIsRHSByTowerNo(noEta); + paramBase->SetCurrentTowerNum(noEta); + paramBase->init(); + + return paramBase; +} + +} +} diff --git a/detectorSegmentations/src/GridEta_k4geo.cpp b/detectorSegmentations/src/GridEta_k4geo.cpp index fada906e2..78573447e 100644 --- a/detectorSegmentations/src/GridEta_k4geo.cpp +++ b/detectorSegmentations/src/GridEta_k4geo.cpp @@ -1,5 +1,4 @@ #include "detectorSegmentations/GridEta_k4geo.h" -#include "DD4hep/Factories.h" namespace dd4hep { namespace DDSegmentation { diff --git a/detectorSegmentations/src/GridSimplifiedDriftChamber_k4geo.cpp b/detectorSegmentations/src/GridSimplifiedDriftChamber_k4geo.cpp index c13dbbddd..635c1d219 100644 --- a/detectorSegmentations/src/GridSimplifiedDriftChamber_k4geo.cpp +++ b/detectorSegmentations/src/GridSimplifiedDriftChamber_k4geo.cpp @@ -1,5 +1,4 @@ #include "detectorSegmentations/GridSimplifiedDriftChamber_k4geo.h" -#include "DD4hep/Factories.h" namespace dd4hep { namespace DDSegmentation { diff --git a/detectorSegmentations/src/plugins/SegmentationFactories.cpp b/detectorSegmentations/src/plugins/SegmentationFactories.cpp index 4dd15369a..d362660de 100644 --- a/detectorSegmentations/src/plugins/SegmentationFactories.cpp +++ b/detectorSegmentations/src/plugins/SegmentationFactories.cpp @@ -29,4 +29,8 @@ DECLARE_SEGMENTATION(GridRPhiEta_k4geo, create_segmentation) +#include "detectorSegmentations/GridDRcalo_k4geo.h" +DECLARE_SEGMENTATION(GridDRcalo_k4geo, create_segmentation) +#include "detectorSegmentations/FCCSWEndcapTurbine_k4geo.h" +DECLARE_SEGMENTATION(FCCSWEndcapTurbine_k4geo, create_segmentation) diff --git a/example/SteeringFile_IDEA_o1_v03.py b/example/SteeringFile_IDEA_o1_v03.py new file mode 100644 index 000000000..f9adeb07e --- /dev/null +++ b/example/SteeringFile_IDEA_o1_v03.py @@ -0,0 +1,651 @@ +from DDSim.DD4hepSimulation import DD4hepSimulation +from g4units import mm, GeV, MeV +SIM = DD4hepSimulation() + +## The compact XML file, or multiple compact files, if the last one is the closer. +SIM.compactFile = ['../FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml'] +## Lorentz boost for the crossing angle, in radian! +SIM.crossingAngleBoost = 0.0 +SIM.enableDetailedShowerMode = False +SIM.enableG4GPS = False +SIM.enableG4Gun = False +SIM.enableGun = True +## InputFiles for simulation .stdhep, .slcio, .HEPEvt, .hepevt, .pairs, .hepmc, .hepmc.gz, .hepmc.xz, .hepmc.bz2, .hepmc3, .hepmc3.gz, .hepmc3.xz, .hepmc3.bz2, .hepmc3.tree.root files are supported +SIM.inputFiles = [] +## Macro file to execute for runType 'run' or 'vis' +SIM.macroFile = "" +## number of events to simulate, used in batch mode +SIM.numberOfEvents = 1 +## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported +SIM.outputFile = "testIDEA_o1_v03.root" +## Physics list to use in simulation +SIM.physicsList = None +## Verbosity use integers from 1(most) to 7(least) verbose +## or strings: VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL, ALWAYS +SIM.printLevel = 3 +## The type of action to do in this invocation +## batch: just simulate some events, needs numberOfEvents, and input file or gun +## vis: enable visualisation, run the macroFile if it is set +## qt: enable visualisation in Qt shell, run the macroFile if it is set +## run: run the macroFile and exit +## shell: enable interactive session +SIM.runType = "batch" +## Skip first N events when reading a file +SIM.skipNEvents = 0 +## Steering file to change default behaviour +SIM.steeringFile = None +## FourVector of translation for the Smearing of the Vertex position: x y z t +SIM.vertexOffset = [0.0, 0.0, 0.0, 0.0] +## FourVector of the Sigma for the Smearing of the Vertex position: x y z t +SIM.vertexSigma = [0.0, 0.0, 0.0, 0.0] + + +################################################################################ +## Helper holding sensitive detector and other actions. +## +## The default tracker and calorimeter sensitive actions can be set with +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.tracker=('Geant4TrackerWeightedAction', {'HitPositionCombination': 2, 'CollectSingleDeposits': False}) +## >>> SIM.action.calo = "Geant4CalorimeterAction" +## +## The default sensitive actions for calorimeters and trackers are applied based on the sensitive type. +## The list of sensitive types can be changed with +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.trackerSDTypes = ['tracker', 'myTrackerSensType'] +## >>> SIM.calor.calorimeterSDTypes = ['calorimeter', 'myCaloSensType'] +## +## For specific subdetectors specific sensitive detectors can be set based on patterns in the name of the subdetector. +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.mapActions['tpc'] = "TPCSDAction" +## +## and additional parameters for the sensitive detectors can be set when the map is given a tuple +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.mapActions['ecal'] =( "CaloPreShowerSDAction", {"FirstLayerNumber": 1} ) +## +## Additional actions can be set as well with the following syntax variations: +## +## >>> SIM = DD4hepSimulation() +## # single action by name only: +## >>> SIM.action.run = "Geant4TestRunAction" +## # multiple actions with comma-separated names: +## >>> SIM.action.event = "Geant4TestEventAction/Action0,Geant4TestEventAction/Action1" +## # single action by tuple of name and parameter dict: +## >>> SIM.action.track = ( "Geant4TestTrackAction", {"Property_int": 10} ) +## # single action by dict of name and parameter dict: +## >>> SIM.action.step = { "name": "Geant4TestStepAction", "parameter": {"Property_int": 10} } +## # multiple actions by list of dict of name and parameter dict: +## >>> SIM.action.stack = [ { "name": "Geant4TestStackAction", "parameter": {"Property_int": 10} } ] +## +## On the command line or in python, these actions can be specified as JSON strings: +## $ ddsim --action.stack '{ "name": "Geant4TestStackAction", "parameter": { "Property_int": 10 } }' +## or +## >>> SIM.action.stack = ''' +## { +## "name": "Geant4TestStackAction", +## "parameter": { +## "Property_int": 10, +## "Property_double": "1.0*mm" +## } +## } +## ''' +## +## +################################################################################ + +## set the default calorimeter action +SIM.action.calo = "Geant4ScintillatorCalorimeterAction" + +## List of patterns matching sensitive detectors of type Calorimeter. +SIM.action.calorimeterSDTypes = ['calorimeter', 'DRcaloSiPMSD'] + +## set the default event action +SIM.action.event = [] + +## Create a map of patterns and actions to be applied to sensitive detectors. +## +## Example: if the name of the detector matches 'tpc' the TPCSDAction is used. +## +## SIM.action.mapActions['tpc'] = "TPCSDAction" +## +SIM.action.mapActions = {'DRcalo': 'DRCaloSDAction'} + +## set the default run action +SIM.action.run = [] + +## set the default stack action +SIM.action.stack = [] + +## set the default step action +SIM.action.step = [] + +## set the default track action +SIM.action.track = [] + +## set the default tracker action +SIM.action.tracker = ('Geant4TrackerWeightedAction', {'HitPositionCombination': 2, 'CollectSingleDeposits': False}) + +## List of patterns matching sensitive detectors of type Tracker. +SIM.action.trackerSDTypes = ['tracker'] + + +################################################################################ +## Configuration for the magnetic field (stepper) +################################################################################ +SIM.field.delta_chord = 0.25 +SIM.field.delta_intersection = 0.001 +SIM.field.delta_one_step = 0.01 +SIM.field.eps_max = 0.001 +SIM.field.eps_min = 5e-05 +SIM.field.equation = "Mag_UsualEqRhs" +SIM.field.largest_step = 10000.0 +SIM.field.min_chord_step = 0.01 +SIM.field.stepper = "ClassicalRK4" + + +################################################################################ +## Configuration for sensitive detector filters +## +## Set the default filter for 'tracker' +## >>> SIM.filter.tracker = "edep1kev" +## Use no filter for 'calorimeter' by default +## >>> SIM.filter.calo = "" +## +## Assign a filter to a sensitive detector via pattern matching +## >>> SIM.filter.mapDetFilter['FTD'] = "edep1kev" +## +## Or more than one filter: +## >>> SIM.filter.mapDetFilter['FTD'] = ["edep1kev", "geantino"] +## +## Don't use the default filter or anything else: +## >>> SIM.filter.mapDetFilter['TPC'] = None ## or "" or [] +## +## Create a custom filter. The dictionary is used to instantiate the filter later on +## >>> SIM.filter.filters['edep3kev'] = dict(name="EnergyDepositMinimumCut/3keV", parameter={"Cut": 3.0*keV} ) +## +## +################################################################################ + +## +## default filter for calorimeter sensitive detectors; +## this is applied if no other filter is used for a calorimeter +## +SIM.filter.calo = "edep0" + +## list of filter objects: map between name and parameter dictionary +SIM.filter.filters = {'geantino': {'name': 'GeantinoRejectFilter/GeantinoRejector', 'parameter': {}}, 'edep1kev': {'name': 'EnergyDepositMinimumCut', 'parameter': {'Cut': 0.001}}, 'edep0': {'name': 'EnergyDepositMinimumCut/Cut0', 'parameter': {'Cut': 0.0}}} + +## a map between patterns and filter objects, using patterns to attach filters to sensitive detector +SIM.filter.mapDetFilter = {} + +## default filter for tracking sensitive detectors; this is applied if no other filter is used for a tracker +SIM.filter.tracker = "edep1kev" + + +################################################################################ +## Configuration for the Detector Construction. +################################################################################ +SIM.geometry.dumpGDML = "" +SIM.geometry.dumpHierarchy = 0 + +## Print Debug information about Elements +SIM.geometry.enableDebugElements = False + +## Print Debug information about Materials +SIM.geometry.enableDebugMaterials = False + +## Print Debug information about Placements +SIM.geometry.enableDebugPlacements = False + +## Print Debug information about Reflections +SIM.geometry.enableDebugReflections = False + +## Print Debug information about Regions +SIM.geometry.enableDebugRegions = False + +## Print Debug information about Shapes +SIM.geometry.enableDebugShapes = False + +## Print Debug information about Surfaces +SIM.geometry.enableDebugSurfaces = False + +## Print Debug information about Volumes +SIM.geometry.enableDebugVolumes = False + +## Print information about placements +SIM.geometry.enablePrintPlacements = False + +## Print information about Sensitives +SIM.geometry.enablePrintSensitives = False + + +################################################################################ +## Configuration for the GuineaPig InputFiles +################################################################################ + +## Set the number of pair particles to simulate per event. +## Only used if inputFile ends with ".pairs" +## If "-1" all particles will be simulated in a single event +## +SIM.guineapig.particlesPerEvent = "-1" + + +################################################################################ +## Configuration for the DDG4 ParticleGun +################################################################################ + +## direction of the particle gun, 3 vector +SIM.gun.direction = (1., 1., 1.) + +## choose the distribution of the random direction for theta +## +## Options for random distributions: +## +## 'uniform' is the default distribution, flat in theta +## 'cos(theta)' is flat in cos(theta) +## 'eta', or 'pseudorapidity' is flat in pseudorapity +## 'ffbar' is distributed according to 1+cos^2(theta) +## +## Setting a distribution will set isotrop = True +## +SIM.gun.distribution = None + +## Total energy (including mass) for the particle gun. +## +## If not None, it will overwrite the setting of momentumMin and momentumMax +SIM.gun.energy = None + +## Maximal pseudorapidity for random distibution (overrides thetaMin) +SIM.gun.etaMax = None + +## Minimal pseudorapidity for random distibution (overrides thetaMax) +SIM.gun.etaMin = None + +## isotropic distribution for the particle gun +## +## use the options phiMin, phiMax, thetaMin, and thetaMax to limit the range of randomly distributed directions +## if one of these options is not None the random distribution will be set to True and cannot be turned off! +## +SIM.gun.isotrop = False + +## Maximal momentum when using distribution (default = 0.0) +SIM.gun.momentumMax = 10000.0 + +## Minimal momentum when using distribution (default = 0.0) +SIM.gun.momentumMin = 0.0 +SIM.gun.multiplicity = 1 +SIM.gun.particle = "mu-" + +## Maximal azimuthal angle for random distribution +SIM.gun.phiMax = None + +## Minimal azimuthal angle for random distribution +SIM.gun.phiMin = None + +## position of the particle gun, 3 vector +SIM.gun.position = (0.0, 0.0, 0.0) + +## Maximal polar angle for random distribution +SIM.gun.thetaMax = None + +## Minimal polar angle for random distribution +SIM.gun.thetaMin = None + + +################################################################################ +## Configuration for the hepmc3 InputFiles +################################################################################ + +## Set the name of the attribute contraining color flow information index 0. +SIM.hepmc3.Flow1 = "flow1" + +## Set the name of the attribute contraining color flow information index 1. +SIM.hepmc3.Flow2 = "flow2" + +## Set to false if the input should be opened with the hepmc2 ascii reader. +## +## If ``True`` a '.hepmc' file will be opened with the HEPMC3 Reader Factory. +## +## Defaults to true if DD4hep was build with HEPMC3 support. +## +SIM.hepmc3.useHepMC3 = True + + +################################################################################ +## Configuration for Input Files. +################################################################################ + +## Set one or more functions to configure input steps. +## +## The functions must take a ``DD4hepSimulation`` object as their only argument and return the created generatorAction +## ``gen`` (for example). +## +## For example one can add this to the ddsim steering file: +## +## def exampleUserPlugin(dd4hepSimulation): +## '''Example code for user created plugin. +## +## :param DD4hepSimulation dd4hepSimulation: The DD4hepSimulation instance, so all parameters can be accessed +## :return: GeneratorAction +## ''' +## from DDG4 import GeneratorAction, Kernel +## # Geant4InputAction is the type of plugin, Cry1 just an identifier +## gen = GeneratorAction(Kernel(), 'Geant4InputAction/Cry1' , True) +## # CRYEventReader is the actual plugin, steeringFile its constructor parameter +## gen.Input = 'CRYEventReader|' + 'steeringFile' +## # we can give a dictionary of Parameters that has to be interpreted by the setParameters function of the plugin +## gen.Parameters = {'DataFilePath': '/path/to/files/data'} +## gen.enableUI() +## return gen +## +## SIM.inputConfig.userInputPlugin = exampleUserPlugin +## +## Repeat function definition and assignment to add multiple input steps +## +## +SIM.inputConfig.userInputPlugin = [] + + +################################################################################ +## Configuration for the generator-level InputFiles +################################################################################ + +## Set the name of the collection containing the MCParticle input. +## Default is "MCParticle". +## +SIM.lcio.mcParticleCollectionName = "MCParticle" + + +################################################################################ +## Configuration for the LCIO output file settings +################################################################################ + +## The event number offset to write in slcio output file. E.g setting it to 42 will start counting events from 42 instead of 0 +SIM.meta.eventNumberOffset = 0 + +## Event parameters to write in every event. Use C/F/I ids to specify parameter type. E.g parameterName/F=0.42 to set a float parameter +SIM.meta.eventParameters = [] + +## The run number offset to write in slcio output file. E.g setting it to 42 will start counting runs from 42 instead of 0 +SIM.meta.runNumberOffset = 0 + + +################################################################################ +## Configuration for the output levels of DDG4 components +################################################################################ + +## Output level for geometry. +SIM.output.geometry = 2 + +## Output level for input sources +SIM.output.inputStage = 3 + +## Output level for Geant4 kernel +SIM.output.kernel = 3 + +## Output level for ParticleHandler +SIM.output.part = 3 + +## Output level for Random Number Generator setup +SIM.output.random = 6 + + +################################################################################ +## Configuration for Output Files. +################################################################################ + +## Use the DD4HEP output plugin regardless of outputfilename. +SIM.outputConfig.forceDD4HEP = False + +## Use the EDM4HEP output plugin regardless of outputfilename. +SIM.outputConfig.forceEDM4HEP = False + +## Use the LCIO output plugin regardless of outputfilename. +SIM.outputConfig.forceLCIO = False + +## Set a function to configure the outputFile. +## +## The function must take a ``DD4hepSimulation`` object as its only argument and return ``None``. +## +## For example one can add this to the ddsim steering file: +## +## def exampleUserPlugin(dd4hepSimulation): +## '''Example code for user created plugin. +## +## :param DD4hepSimulation dd4hepSimulation: The DD4hepSimulation instance, so all parameters can be accessed +## :return: None +## ''' +## from DDG4 import EventAction, Kernel +## dd = dd4hepSimulation # just shorter variable name +## evt_root = EventAction(Kernel(), 'Geant4Output2ROOT/' + dd.outputFile, True) +## evt_root.HandleMCTruth = True or False +## evt_root.Control = True +## output = dd.outputFile +## if not dd.outputFile.endswith(dd.outputConfig.myExtension): +## output = dd.outputFile + dd.outputConfig.myExtension +## evt_root.Output = output +## evt_root.enableUI() +## Kernel().eventAction().add(evt_root) +## return None +## +## SIM.outputConfig.userOutputPlugin = exampleUserPlugin +## # arbitrary options can be created and set via the steering file or command line +## SIM.outputConfig.myExtension = '.csv' +## + +def Geant4Output2EDM4hep_DRC_plugin(dd4hepSimulation): + from DDG4 import EventAction, Kernel + evt_root = EventAction(Kernel(), 'Geant4Output2EDM4hep_DRC/' + dd4hepSimulation.outputFile, True) + evt_root.Control = True + output = dd4hepSimulation.outputFile + evt_root.Output = output + evt_root.enableUI() + Kernel().eventAction().add(evt_root) + return None + +SIM.outputConfig.userOutputPlugin = Geant4Output2EDM4hep_DRC_plugin + +################################################################################ +## Configuration for the Particle Handler/ MCTruth treatment +################################################################################ + +## Enable lots of printout on simulated hits and MC-truth information +SIM.part.enableDetailedHitsAndParticleInfo = False + +## Keep all created particles +SIM.part.keepAllParticles = False + +## Minimal distance between particle vertex and endpoint of parent after +## which the vertexIsNotEndpointOfParent flag is set +## +SIM.part.minDistToParentVertex = 2.2e-14 + +## MinimalKineticEnergy to store particles created in the tracking region +SIM.part.minimalKineticEnergy = 1.0 + +## Printout at End of Tracking +SIM.part.printEndTracking = False + +## Printout at Start of Tracking +SIM.part.printStartTracking = False + +## List of processes to save, on command line give as whitespace separated string in quotation marks +SIM.part.saveProcesses = ['Decay'] + +## Optionally enable an extended Particle Handler +SIM.part.userParticleHandler = "Geant4TCUserParticleHandler" + + +################################################################################ +## Configuration for the PhysicsList and Monte Carlo particle selection. +## +## To load arbitrary plugins, add a function to be executed. +## +## The function must take the DDG4.Kernel() object as the only argument. +## +## For example, add a function definition and the call to a steering file:: +## +## def setupCerenkov(kernel): +## from DDG4 import PhysicsList +## seq = kernel.physicsList() +## cerenkov = PhysicsList(kernel, 'Geant4CerenkovPhysics/CerenkovPhys') +## cerenkov.MaxNumPhotonsPerStep = 10 +## cerenkov.MaxBetaChangePerStep = 10.0 +## cerenkov.TrackSecondariesFirst = True +## cerenkov.VerboseLevel = 2 +## cerenkov.enableUI() +## seq.adopt(cerenkov) +## ph = PhysicsList(kernel, 'Geant4OpticalPhotonPhysics/OpticalGammaPhys') +## ph.addParticleConstructor('G4OpticalPhoton') +## ph.VerboseLevel = 2 +## ph.enableUI() +## seq.adopt(ph) +## return None +## +## SIM.physics.setupUserPhysics(setupCerenkov) +## +## # End of example +## +################################################################################ + +## Set of Generator Statuses that are used to mark unstable particles that should decay inside of Geant4. +## +SIM.physics.alternativeDecayStatuses = set() + +## If true, add decay processes for all particles. +## +## Only enable when creating a physics list not based on an existing Geant4 list! +## +SIM.physics.decays = False + +## The name of the Geant4 Physics list. +SIM.physics.list = "FTFP_BERT" + +## location of particle.tbl file containing extra particles and their lifetime information +## +## For example in $DD4HEP/examples/DDG4/examples/particle.tbl +## +SIM.physics.pdgfile = None + +## The global geant4 rangecut for secondary production +## +## Default is 0.7 mm as is the case in geant4 10 +## +## To disable this plugin and be absolutely sure to use the Geant4 default range cut use "None" +## +## Set printlevel to DEBUG to see a printout of all range cuts, +## but this only works if range cut is not "None" +## +SIM.physics.rangecut = 0.7 + +## Set of PDG IDs that will not be passed from the input record to Geant4. +## +## Quarks, gluons and W's Z's etc should not be treated by Geant4 +## +SIM.physics.rejectPDGs = {1, 2, 3, 4, 5, 6, 3201, 3203, 4101, 4103, 21, 23, 24, 5401, 25, 2203, 5403, 3101, 3103, 4403, 2101, 5301, 2103, 5303, 4301, 1103, 4303, 5201, 5203, 3303, 4201, 4203, 5101, 5103, 5503} + +## Set of PDG IDs for particles that should not be passed to Geant4 if their properTime is 0. +## +## The properTime of 0 indicates a documentation to add FSR to a lepton for example. +## +SIM.physics.zeroTimePDGs = {17, 11, 13, 15} + + +def setupOpticalPhysics(kernel): + from DDG4 import PhysicsList + seq = kernel.physicsList() + cerenkov = PhysicsList(kernel, 'Geant4CerenkovPhysics/CerenkovPhys') + cerenkov.TrackSecondariesFirst = True + cerenkov.VerboseLevel = 1 + cerenkov.enableUI() + seq.adopt(cerenkov) + + scint = PhysicsList(kernel, 'Geant4ScintillationPhysics/ScintillationPhys') + scint.VerboseLevel = 1 + scint.TrackSecondariesFirst = True + scint.BoundaryInvokeSD = True + scint.enableUI() + seq.adopt(scint) + + opt = PhysicsList(kernel, 'Geant4OpticalPhotonPhysics/OpticalGammaPhys') + opt.addParticleConstructor('G4OpticalPhoton') + opt.VerboseLevel = 1 + opt.BoundaryInvokeSD = True + opt.enableUI() + seq.adopt(opt) + + return None +SIM.physics.setupUserPhysics(setupOpticalPhysics) + +def setupDRCFastSim(kernel): + from DDG4 import DetectorConstruction, Geant4, PhysicsList + geant4 = Geant4(kernel) + seq = geant4.detectorConstruction() + # Create a model for fast simulation + model = DetectorConstruction(kernel, str("Geant4DRCFiberModel/ShowerModel") ) + # Mandatory model parameters + model.RegionName = 'FastSimOpFiberRegion' + model.Enable = True + model.ApplicableParticles = ['opticalphoton'] + model.enableUI() + seq.adopt(model) + # Now build the physics list: + phys = kernel.physicsList() + ph = PhysicsList(kernel, str('Geant4FastPhysics/FastPhysicsList')) + ph.EnabledParticles = ['opticalphoton'] + ph.BeVerbose = True + ph.enableUI() + phys.adopt(ph) + phys.dump() +SIM.physics.setupUserPhysics(setupDRCFastSim) + +################################################################################ +## Properties for the random number generator +################################################################################ + +## If True, calculate random seed for each event basedon eventID and runID +## Allows reproducibility even whenSkippingEvents +SIM.random.enableEventSeed = False +SIM.random.file = None +SIM.random.luxury = 1 +SIM.random.replace_gRandom = True +SIM.random.seed = None +SIM.random.type = None + + +################################################################################ +## Configuration for setting commands to run during different phases. +## +## In this section, one can configure commands that should be run during the different phases of the Geant4 execution. +## +## 1. Configuration +## 2. Initialization +## 3. Pre Run +## 4. Post Run +## 5. Terminate / Finalization +## +## For example, one can add +## +## >>> SIM.ui.commandsConfigure = ['/physics_lists/em/SyncRadiation true'] +## +## Further details should be taken from the Geant4 documentation. +## +################################################################################ + +## List of UI commands to run during the 'Configure' phase. +SIM.ui.commandsConfigure = [] + +## List of UI commands to run during the 'Initialize' phase. +SIM.ui.commandsInitialize = [] + +## List of UI commands to run during the 'PostRun' phase. +SIM.ui.commandsPostRun = [] + +## List of UI commands to run during the 'PreRun' phase. +SIM.ui.commandsPreRun = [] + +## List of UI commands to run during the 'Terminate' phase. +SIM.ui.commandsTerminate = [] diff --git a/plugins/CaloPreShowerSDAction.cpp b/plugins/CaloPreShowerSDAction.cpp index fdefaf205..4f8bcee38 100644 --- a/plugins/CaloPreShowerSDAction.cpp +++ b/plugins/CaloPreShowerSDAction.cpp @@ -1,9 +1,5 @@ #include "DD4hep/Version.h" #include "DDG4/Geant4SensDetAction.inl" -#include "DDG4/Geant4EventAction.h" -#include "DDG4/Geant4Mapping.h" -#include "G4OpticalPhoton.hh" -#include "G4VProcess.hh" #if DD4HEP_VERSION_GE(1, 21) #define GEANT4_CONST_STEP const @@ -132,4 +128,5 @@ namespace dd4hep { #include "DDG4/Factories.h" + DECLARE_GEANT4SENSITIVE( CaloPreShowerSDAction ) diff --git a/plugins/DRCaloFastSimModel.cpp b/plugins/DRCaloFastSimModel.cpp new file mode 100644 index 000000000..5caa21456 --- /dev/null +++ b/plugins/DRCaloFastSimModel.cpp @@ -0,0 +1,376 @@ +// Framework include files +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// C/C++ include files +#include "DRCaloFastSimModel.h" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep +{ + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim + { + + class DRCFiberModel + { + public: + // G4FastSimHitMaker hitMaker { }; + G4OpBoundaryProcess *pOpBoundaryProc{nullptr}; + G4OpAbsorption *pOpAbsorption{nullptr}; + G4OpWLS *pOpWLS{nullptr}; + G4bool fProcAssigned{false}; + + FastFiberData mDataPrevious{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))}; + FastFiberData mDataCurrent{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))}; + + G4int fSafety{2}; + G4double mNtransport{0.}; + G4double mTransportUnit{0.}; + G4ThreeVector mFiberPos{G4ThreeVector(0)}; + G4ThreeVector mFiberAxis{G4ThreeVector(0)}; + G4bool fKill{false}; + G4bool fTransported{false}; + G4bool fSwitch{true}; + G4int fVerbose{0}; + + G4bool checkTotalInternalReflection(const G4Track *track) + { + if (!fProcAssigned) + setPostStepProc(track); // locate OpBoundaryProcess only once + + if (track->GetTrackStatus() == fStopButAlive || track->GetTrackStatus() == fStopAndKill) + return false; + + // accumulate step length + mDataCurrent.AddStepLengthInterval(track->GetStepLength()); + + G4int theStatus = pOpBoundaryProc->GetStatus(); + + if (fVerbose > 1) + { + G4cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID(); + G4cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus; + G4cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << G4endl; + } + + // skip exceptional iteration with FresnelReflection + if (theStatus == G4OpBoundaryProcessStatus::FresnelReflection) + mDataCurrent.SetOpBoundaryStatus(theStatus); + + // some cases have a status StepTooSmall when the reflection happens between the boundary of cladding & outer volume (outside->cladding) since the outer volume is not a G4Region + if (theStatus == G4OpBoundaryProcessStatus::TotalInternalReflection || theStatus == G4OpBoundaryProcessStatus::StepTooSmall) + { + if (theStatus != G4OpBoundaryProcessStatus::TotalInternalReflection) + { // skip StepTooSmall if the track already has TotalInternalReflection + if (mDataCurrent.GetOpBoundaryStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection) + return false; + if (mDataPrevious.GetOpBoundaryStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection) + return false; + } + + G4int trackID = track->GetTrackID(); + G4double kineticEnergy = track->GetKineticEnergy(); + G4double globalTime = track->GetGlobalTime(); + G4double pathLength = track->GetStepLength(); + G4ThreeVector globalPosition = track->GetPosition(); + G4ThreeVector momentumDirection = track->GetMomentumDirection(); + G4ThreeVector polarization = track->GetPolarization(); + + auto candidate = FastFiberData(trackID, kineticEnergy, globalTime, pathLength, globalPosition, momentumDirection, polarization, theStatus); + if (pOpAbsorption != nullptr) + candidate.SetAbsorptionNILL(pOpAbsorption->GetNumberOfInteractionLengthLeft()); + if (pOpWLS != nullptr) + candidate.SetWLSNILL(pOpWLS->GetNumberOfInteractionLengthLeft()); + + G4bool repetitive = false; + if (candidate.checkRepetitive(mDataCurrent, false) && mDataCurrent.checkRepetitive(mDataPrevious)) + repetitive = true; + + mDataPrevious = mDataCurrent; + mDataCurrent = candidate; + + return repetitive; + } + + return false; + } + + G4bool checkAbsorption(const G4double prevNILL, const G4double currentNILL) + { + if (prevNILL < 0. || currentNILL < 0.) + return false; // the number of interaction length left has to be reset + if (prevNILL == currentNILL) + return false; // no absorption + if (prevNILL == DBL_MAX || currentNILL == DBL_MAX) + return false; // NILL is re-initialized + + G4double deltaNILL = prevNILL - currentNILL; + + if (currentNILL - deltaNILL * (mNtransport + fSafety) < 0.) + return true; // absorbed before reaching fiber end + + return false; + } + + G4bool checkNILL() + { + if (!fTransported) + return true; // do nothing if the track is not already transported + + G4double wlsNILL = DBL_MAX; + G4double absorptionNILL = DBL_MAX; + + if (pOpWLS != nullptr) + { + wlsNILL = pOpWLS->GetNumberOfInteractionLengthLeft(); + if (mDataPrevious.GetWLSNILL() == DBL_MAX || mDataCurrent.GetWLSNILL() == DBL_MAX) + return true; // NILL is re-initialized + } + + if (pOpAbsorption != nullptr) + { + absorptionNILL = pOpAbsorption->GetNumberOfInteractionLengthLeft(); + if (mDataPrevious.GetAbsorptionNILL() == DBL_MAX || mDataCurrent.GetAbsorptionNILL() == DBL_MAX) + return true; // NILL is re-initialized + } + + if (wlsNILL < 0. || absorptionNILL < 0.) + return true; // let GEANT4 to reset them + + G4double deltaWlsNILL = mDataPrevious.GetWLSNILL() - mDataCurrent.GetWLSNILL(); + G4double deltaAbsorptionNILL = mDataPrevious.GetAbsorptionNILL() - mDataCurrent.GetAbsorptionNILL(); + + G4double finalWlsNILL = wlsNILL - deltaWlsNILL * fSafety; + G4double finalAbsorptionNILL = absorptionNILL - deltaAbsorptionNILL * fSafety; + + // prevent double counting of the probability of getting absorbed (which already estimated before transportation) + // reset NILL again + if (finalWlsNILL < 0. || finalAbsorptionNILL < 0.) + return false; + + return true; + } + + void setPostStepProc(const G4Track *track) + { + G4ProcessManager *pm = track->GetDefinition()->GetProcessManager(); + auto postStepProcessVector = pm->GetPostStepProcessVector(); + + for (unsigned int np = 0; np < postStepProcessVector->entries(); np++) + { + auto theProcess = (*postStepProcessVector)[np]; + + auto theType = theProcess->GetProcessType(); + + if (theType != fOptical) + continue; + + if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpBoundary) + pOpBoundaryProc = dynamic_cast(theProcess); + else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpAbsorption) + pOpAbsorption = dynamic_cast(theProcess); + else if (theProcess->GetProcessSubType() == G4OpProcessSubType::fOpWLS) + pOpWLS = dynamic_cast(theProcess); + } + + fProcAssigned = true; + + return; + } + + void reset() + { + mNtransport = 0.; + mTransportUnit = 0.; + mFiberPos = G4ThreeVector(0); + mFiberAxis = G4ThreeVector(0); + fKill = false; + fTransported = false; + mDataCurrent.reset(); + mDataPrevious.reset(); + } + + void print() + { + if (fVerbose > 1) + { + G4cout << G4endl; + + G4cout << "mDataPrevious.trackID = " << mDataPrevious.trackID; + G4cout << " | .mOpBoundaryStatus = " << std::setw(4) << mDataPrevious.GetOpBoundaryStatus(); + G4cout << " | .mStepLengthInterval = " << mDataPrevious.GetStepLengthInterval() << G4endl; + + if (fVerbose > 2) + { + G4cout << " | globalPosition = (" << std::setw(9) << mDataPrevious.globalPosition.x(); + G4cout << "," << std::setw(9) << mDataPrevious.globalPosition.y(); + G4cout << "," << std::setw(9) << mDataPrevious.globalPosition.z() << ")" << G4endl; + + G4cout << " | momentumDirection = (" << std::setw(9) << mDataPrevious.momentumDirection.x(); + G4cout << "," << std::setw(9) << mDataPrevious.momentumDirection.y(); + G4cout << "," << std::setw(9) << mDataPrevious.momentumDirection.z() << ")" << G4endl; + + G4cout << " | polarization = (" << std::setw(9) << mDataPrevious.polarization.x(); + G4cout << "," << std::setw(9) << mDataPrevious.polarization.y(); + G4cout << "," << std::setw(9) << mDataPrevious.polarization.z() << ")" << G4endl; + + G4cout << " | globalTime = " << std::setw(9) << mDataPrevious.globalTime << G4endl; + G4cout << " | WLSNILL = " << std::setw(9) << mDataPrevious.GetWLSNILL() << G4endl; + G4cout << " | AbsorptionNILL = " << std::setw(9) << mDataPrevious.GetAbsorptionNILL() << G4endl; + } + + G4cout << "mDataCurrent.trackID = " << mDataCurrent.trackID; + G4cout << " | .mOpBoundaryStatus = " << std::setw(4) << mDataCurrent.GetOpBoundaryStatus() << G4endl; + + if (fVerbose > 2) + { + G4cout << " | globalPosition = (" << std::setw(9) << mDataCurrent.globalPosition.x(); + G4cout << "," << std::setw(9) << mDataCurrent.globalPosition.y(); + G4cout << "," << std::setw(9) << mDataCurrent.globalPosition.z() << ")" << G4endl; + + G4cout << " | momentumDirection = (" << std::setw(9) << mDataCurrent.momentumDirection.x(); + G4cout << "," << std::setw(9) << mDataCurrent.momentumDirection.y(); + G4cout << "," << std::setw(9) << mDataCurrent.momentumDirection.z() << ")" << G4endl; + + G4cout << " | polarization = (" << std::setw(9) << mDataCurrent.polarization.x(); + G4cout << "," << std::setw(9) << mDataCurrent.polarization.y(); + G4cout << "," << std::setw(9) << mDataCurrent.polarization.z() << ")" << G4endl; + + G4cout << " | globalTime = " << std::setw(9) << mDataCurrent.globalTime << G4endl; + G4cout << " | WLSNILL = " << std::setw(9) << mDataCurrent.GetWLSNILL() << G4endl; + G4cout << " | AbsorptionNILL = " << std::setw(9) << mDataCurrent.GetAbsorptionNILL() << G4endl; + } + + G4cout << G4endl; + } + } + }; + + template <> + void Geant4FSShowerModel::initialize() + { + this->m_applicablePartNames.emplace_back("opticalphoton"); + } + + template <> + void Geant4FSShowerModel::constructSensitives(Geant4DetectorConstructionContext *ctxt) + { + this->Geant4FastSimShowerModel::constructSensitives(ctxt); + } + + template <> + void Geant4FSShowerModel::modelShower(const G4FastTrack &fasttrack, G4FastStep &faststep) + { + auto *track = fasttrack.GetPrimaryTrack(); + + if (locals.fKill) + { // absorption + faststep.ProposeTotalEnergyDeposited(track->GetKineticEnergy()); + faststep.KillPrimaryTrack(); + + return; + } + + if (locals.fTransported) + return; // reset NILL if the track did not meet NILL check + + double timeUnit = locals.mDataCurrent.globalTime - locals.mDataPrevious.globalTime; + auto posShift = locals.mTransportUnit * locals.mNtransport * locals.mFiberAxis; // #TODO apply shift for xy direction as well + double timeShift = timeUnit * locals.mNtransport; + + faststep.ProposePrimaryTrackFinalPosition(track->GetPosition() + posShift, false); + faststep.ProposePrimaryTrackFinalTime(track->GetGlobalTime() + timeShift); + faststep.ProposePrimaryTrackFinalKineticEnergy(track->GetKineticEnergy()); + faststep.ProposePrimaryTrackFinalMomentumDirection(track->GetMomentumDirection(), false); + faststep.ProposePrimaryTrackFinalPolarization(track->GetPolarization(), false); + locals.fTransported = true; + return; + } + + template <> + bool Geant4FSShowerModel::check_applicability(const G4ParticleDefinition &particle) + { + return &particle == G4OpticalPhoton::OpticalPhotonDefinition(); + } + + template <> + bool Geant4FSShowerModel::check_trigger(const G4FastTrack &fasttrack) + { + if (!locals.fSwitch) + return false; // turn on/off the model + + const G4Track *track = fasttrack.GetPrimaryTrack(); + + // reset when moving to the next track + if (locals.mDataCurrent.trackID != track->GetTrackID()) + locals.reset(); + + // make sure that the track does not get absorbed after transportation, as number of interaction length left is reset when doing transportation + if (!locals.checkNILL()) + return true; // track is already transported but did not pass NILL check, attempt to reset NILL + + if (locals.fTransported) + { // track is already transported and did pass NILL check, nothing to do + if (locals.mFiberAxis.dot(track->GetMomentumDirection()) * locals.mFiberAxis.dot(locals.mDataCurrent.momentumDirection) < 0) // different propagation direction (e.g. mirror) + locals.reset(); + + return false; + } + + if (!locals.checkTotalInternalReflection(track)) + return false; // nothing to do if the track has no repetitive total internal reflection + + auto theTouchable = track->GetTouchableHandle(); + auto solid = theTouchable->GetSolid(); + + if (solid->GetEntityType() != "G4Tubs") + return false; // only works for G4Tubs at the moment + + if (locals.fVerbose > 0) + locals.print(); // at this point, the track should have passed all prerequisites before entering computationally heavy operations + + G4Tubs *tubs = static_cast(solid); + G4double fiberLen = 2. * tubs->GetZHalfLength(); + + locals.mFiberPos = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(G4ThreeVector(0., 0., 0.)); + locals.mFiberAxis = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformAxis(G4ThreeVector(0., 0., 1.)); + + auto delta = locals.mDataCurrent.globalPosition - locals.mDataPrevious.globalPosition; + locals.mTransportUnit = delta.dot(locals.mFiberAxis); + + // estimate the number of expected total internal reflections before reaching fiber end + auto fiberEnd = (locals.mTransportUnit > 0.) ? locals.mFiberPos + locals.mFiberAxis * fiberLen / 2. : locals.mFiberPos - locals.mFiberAxis * fiberLen / 2.; + auto toEnd = fiberEnd - track->GetPosition(); + G4double toEndAxis = toEnd.dot(locals.mFiberAxis); + G4double maxTransport = std::floor(toEndAxis / locals.mTransportUnit); + locals.mNtransport = maxTransport - locals.fSafety; + + if (locals.mNtransport < 1.) + return false; // require at least n = fSafety of total internal reflections at the end + + if (locals.checkAbsorption(locals.mDataPrevious.GetWLSNILL(), locals.mDataCurrent.GetWLSNILL())) + return false; // do nothing if WLS happens before reaching fiber end + if (locals.checkAbsorption(locals.mDataPrevious.GetAbsorptionNILL(), locals.mDataCurrent.GetAbsorptionNILL())) + locals.fKill = true; // absorbed before reaching fiber end + + return true; + } + + typedef Geant4FSShowerModel Geant4DRCFiberModel; + } +} + +#include + +DECLARE_GEANT4ACTION_NS(dd4hep::sim, Geant4DRCFiberModel) diff --git a/plugins/DRCaloFastSimModel.h b/plugins/DRCaloFastSimModel.h new file mode 100644 index 000000000..a8f29cb14 --- /dev/null +++ b/plugins/DRCaloFastSimModel.h @@ -0,0 +1,78 @@ +#include "G4OpBoundaryProcess.hh" + + +struct FastFiberData +{ +public: + FastFiberData(G4int, G4double, G4double, G4double, G4ThreeVector, G4ThreeVector, G4ThreeVector, G4int status = G4OpBoundaryProcessStatus::Undefined); + FastFiberData& operator=(const FastFiberData& theData) = default; + + FastFiberData(const FastFiberData& theData) = default; + ~FastFiberData() = default; + + + void reset(); + + G4double GetAbsorptionNILL() { return mOpAbsorptionNumIntLenLeft; } + void SetAbsorptionNILL(G4double in) { mOpAbsorptionNumIntLenLeft = in; } + + G4double GetWLSNILL() { return mOpWLSNumIntLenLeft; } + void SetWLSNILL(G4double in) { mOpWLSNumIntLenLeft = in; } + + G4int GetOpBoundaryStatus() { return mOpBoundaryStatus; } + void SetOpBoundaryStatus(G4int in) { mOpBoundaryStatus = in; } + + G4double GetStepLengthInterval() { return mStepLengthInterval; } + void AddStepLengthInterval(G4double in) { mStepLengthInterval += in; } + + G4bool checkRepetitive(const FastFiberData, G4bool checkInterval = true); + + G4int trackID; + G4double kineticEnergy; + G4double globalTime; + G4double pathLength; + G4ThreeVector globalPosition; + G4ThreeVector momentumDirection; + G4ThreeVector polarization; + +private: + G4int mOpBoundaryStatus; + G4double mOpAbsorptionNumIntLenLeft; + G4double mOpWLSNumIntLenLeft; + G4double mStepLengthInterval; +}; + +FastFiberData::FastFiberData(G4int id, G4double en, G4double globTime, G4double path, G4ThreeVector pos, G4ThreeVector mom, G4ThreeVector pol, G4int status) +{ + trackID = id; + kineticEnergy = en; + globalTime = globTime; + pathLength = path; + globalPosition = pos; + momentumDirection = mom; + polarization = pol; + mOpBoundaryStatus = status; + mOpAbsorptionNumIntLenLeft = DBL_MAX; + mOpWLSNumIntLenLeft = DBL_MAX; + mStepLengthInterval = 0.; +} + +G4bool FastFiberData::checkRepetitive(const FastFiberData theData, G4bool checkInterval) +{ + if (this->trackID != theData.trackID) + return false; + if (this->mOpBoundaryStatus != theData.mOpBoundaryStatus) + return false; + if (checkInterval && std::abs(this->mStepLengthInterval - theData.mStepLengthInterval) > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) + return false; + + return true; +} + +void FastFiberData::reset() +{ + this->mOpBoundaryStatus = G4OpBoundaryProcessStatus::Undefined; + this->mOpAbsorptionNumIntLenLeft = DBL_MAX; + this->mOpWLSNumIntLenLeft = DBL_MAX; + this->mStepLengthInterval = 0.; +} diff --git a/plugins/FiberDRCaloSDAction.cpp b/plugins/FiberDRCaloSDAction.cpp new file mode 100644 index 000000000..1d946ff21 --- /dev/null +++ b/plugins/FiberDRCaloSDAction.cpp @@ -0,0 +1,257 @@ +// DD4hep Framework include files +#include "DD4hep/Version.h" +#include "DD4hep/Objects.h" +#include "DD4hep/Segmentations.h" +#include "DDG4/Geant4Random.h" +#include "DDG4/Geant4SensDetAction.inl" +#include "DDG4/Geant4Mapping.h" + +#include "G4OpticalPhoton.hh" +// k4geo Framework include files + + +#include "FiberDRCaloSDAction.h" + +#if DD4HEP_VERSION_GE(1, 21) +#define GEANT4_CONST_STEP const +#else +#define GEANT4_CONST_STEP +#endif + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep +{ + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim + { + + /* + * Geant4SensitiveAction sensitive detector for the Dual-readout calorimeter + * + * \author Sungwon Kim + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + + struct DRCData + { + Geant4HitCollection *fHitCollection; + G4int fHCID; + + G4int fWavBin; + G4int fTimeBin; + G4float fWavlenStart; + G4float fWavlenEnd; + G4float fTimeStart; + G4float fTimeEnd; + G4float fWavlenStep; + G4float fTimeStep; + static const int fArraySize = 24; + double fGraph_X[fArraySize] ={1.37760 * dd4hep::eV, + 1.45864 * dd4hep::eV, + 1.54980 * dd4hep::eV, + 1.65312 * dd4hep::eV, + 1.71013 * dd4hep::eV, + 1.77120 * dd4hep::eV, + 1.83680 * dd4hep::eV, + 1.90745 * dd4hep::eV, + 1.98375 * dd4hep::eV, + 2.06640 * dd4hep::eV, + 2.10143 * dd4hep::eV, + 2.13766 * dd4hep::eV, + 2.17516 * dd4hep::eV, + 2.21400 * dd4hep::eV, + 2.25426 * dd4hep::eV, + 2.29600 * dd4hep::eV, + 2.33932 * dd4hep::eV, + 2.38431 * dd4hep::eV, + 2.43106 * dd4hep::eV, + 2.47968 * dd4hep::eV, + 2.53029 * dd4hep::eV, + 2.58300 * dd4hep::eV, + 2.63796 * dd4hep::eV, + 2.69531 * dd4hep::eV}; + double fGraph_Y[fArraySize] ={0.903 * 0.903, + 0.903 * 0.903, + 0.903 * 0.903, + 0.903 * 0.903, + 0.903 * 0.903, + 0.903 * 0.903, + 0.902 * 0.902, + 0.901 * 0.901, + 0.898 * 0.898, + 0.895 * 0.895, + 0.893 * 0.893, + 0.891 * 0.891, + 0.888 * 0.888, + 0.883 * 0.883, + 0.87 * 0.87, + 0.838 * 0.838, + 0.76 * 0.76, + 0.62 * 0.62, + 0.488 * 0.488, + 0.345 * 0.345, + 0.207 * 0.207, + 0.083 * 0.083, + 0.018 * 0.018, + 0.}; + + G4double wavToE(G4double wav) { return CLHEP::h_Planck * CLHEP::c_light / wav; } + + int findWavBin(G4double en) + { + int i = 0; + for (; i < fWavBin + 1; i++) + { + if (en < wavToE((fWavlenStart - static_cast(i) * fWavlenStep) * CLHEP::nm)) + break; + } + + return fWavBin + 1 - i; + } + + int findTimeBin(G4double stepTime) + { + int i = 0; + for (; i < fTimeBin + 1; i++) + { + if (stepTime < ((fTimeStart + static_cast(i) * fTimeStep) * CLHEP::ns)) + break; + } + + return i; + } + + // Linear interpolation function for calculating the efficiency of yellow filter, used in rejectedByYellowFilter + double getFilterEff(G4double G4energy) { + double energy = (G4energy * (dd4hep::MeV / CLHEP::MeV)); // Convert G4 MeV to dd4hep MeV + if (energy <= fGraph_X[0]) // If the photon energy <= 1.37760 eV, than return maximum filter efficiency + return fGraph_Y[0]; + + for(int idx = 0; idx < fArraySize; ++idx) { + if (energy <= fGraph_X[idx]) { + double x1 = fGraph_X[idx - 1]; + double x2 = fGraph_X[idx]; + double y1 = fGraph_Y[idx - 1]; + double y2 = fGraph_Y[idx]; + + return (y1 + ((y2 - y1) / (x2 - x1))*(energy - x1)); // return linear interpolated filter efficiency + } + } + + // This should not happen + std::cout << "Error in Yellow filter efficiency with photon energy : " << energy << " MeV" << std::endl; + std::cout << "Cannot find corresponding filter efficiency" << std::endl; + return 0.; + } + + // If true, then the photon is rejected by yellow filter + bool rejectedByYellowFilter(G4double G4energy, double rndVal) + { + double energy = (G4energy * (dd4hep::MeV / CLHEP::MeV)); // Convert G4 MeV to dd4hep MeV + if ( energy >= fGraph_X[fArraySize-1] ) // Photon w/ E > 2.69531 eV rejected + return true; + + double FilterEff = getFilterEff(G4energy); // Get efficiency of filter using photon's energy + + // filter efficiency == probability of photon accepted by filter + // == Probability of random value (from uniform distribution with range 0 ~ 1) smaller than filter efficiency + // So if the rndVal is larger than the FilterEff, than the photon is rejected + return (rndVal > FilterEff); + } + + DRCData() + : fHitCollection(0), fHCID(-1), fWavBin(120), fTimeBin(650), + fWavlenStart(900.), fWavlenEnd(300.), fTimeStart(5.), fTimeEnd(70.) + { + fWavlenStep = (fWavlenStart - fWavlenEnd) / (float)fWavBin; + fTimeStep = (fTimeEnd - fTimeStart) / (float)fTimeBin; + } + }; + + /// Define collections created by this sensitive action object + template <> + void Geant4SensitiveAction::defineCollections() + { + std::string readout_name = m_sensitive.readout().name(); + m_collectionID = defineCollection(m_sensitive.readout().name()); + std::cout << "defineCollection Geant4DRCalorimeter readout_name : " << readout_name << std::endl; + std::cout << "defineCollection Geant4DRCalorimeter m_collectionID : " << m_collectionID << std::endl; + } + + /// Method for generating hit(s) using the information of G4Step object. + template <> + G4bool Geant4SensitiveAction::process(G4Step GEANT4_CONST_STEP *step, G4TouchableHistory*) + { + if (step->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + return false; + + typedef Geant4DRCalorimeter::Hit Hit; + + Geant4StepHandler h(step); + Geant4HitCollection *coll = collection(m_collectionID); + HitContribution contrib = Hit::extractContribution(step); + + auto theTouchable = step->GetPostStepPoint()->GetTouchable(); + dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager(); + dd4hep::VolumeID volID = volMgr.volumeID(theTouchable); + G4ThreeVector global = step->GetPostStepPoint()->GetPosition(); + G4ThreeVector local = theTouchable->GetHistory()->GetTopTransform().TransformPoint(global); + // MoveUpHistory(GetHistoryDepth - 2) -> tower touchable, local position + dd4hep::Position loc(local.x() * dd4hep::millimeter / CLHEP::millimeter, local.y() * dd4hep::millimeter / CLHEP::millimeter, local.z() * dd4hep::millimeter / CLHEP::millimeter); + dd4hep::Position glob(global.x() * dd4hep::millimeter / CLHEP::millimeter, global.y() * dd4hep::millimeter / CLHEP::millimeter, global.z() * dd4hep::millimeter / CLHEP::millimeter); + + auto cID = m_segmentation->cellID(loc, glob, volID); // This returns cID corresponding to SiPM Wafer + Hit *hit = coll->find(CellIDCompare(cID)); + + G4double hitTime = step->GetPostStepPoint()->GetGlobalTime(); + G4double energy = step->GetTrack()->GetTotalEnergy(); + + // Apply yellow filter here + // Get random number from (0~1) uniform distribution + // If the photon is from Scint. process, calculate the filter eff. based on its energy + // If (random number) > (Eff), the photon is rejected by yellow filter (= do not make hit (or count photon) for that photon) + dd4hep::DDSegmentation::VolumeID ceren = static_cast(m_segmentation->decoder()->get(cID, "c")); + bool IsCeren = static_cast(ceren); + if (!(IsCeren)) + { + Geant4Event& evt = context()->event(); + Geant4Random& rnd = evt.random(); // Get random generator + double rndVal = rnd.uniform(0, 1); // Get random number from uniform distribution [0, 1] + if (m_userData.rejectedByYellowFilter(energy, rndVal)) + return true; + } + + if (!hit) + { + hit = new Geant4DRCalorimeter::Hit(m_userData.fWavlenStep, m_userData.fTimeStep); + hit->cellID = cID; + hit->SetSiPMnum(cID); + hit->SetTimeStart(m_userData.fTimeStart); + hit->SetTimeEnd(m_userData.fTimeEnd); + hit->SetWavlenMax(m_userData.fWavlenStart); + hit->SetWavlenMin(m_userData.fWavlenEnd); + coll->add(hit); + } + + hit->photonCount(); + int wavBin = m_userData.findWavBin(energy); + hit->CountWavlenSpectrum(wavBin); + int timeBin = m_userData.findTimeBin(hitTime); + hit->CountTimeStruct(timeBin); + + hit->position = glob; + hit->energyDeposit += contrib.deposit; + hit->truth.emplace_back(contrib); + + return true; + } + + typedef Geant4SensitiveAction DRCaloSDAction; + } +} + +#include "DDG4/Factories.h" + +DECLARE_GEANT4SENSITIVE(DRCaloSDAction) diff --git a/plugins/FiberDRCaloSDAction.h b/plugins/FiberDRCaloSDAction.h new file mode 100644 index 000000000..11c4e04ad --- /dev/null +++ b/plugins/FiberDRCaloSDAction.h @@ -0,0 +1,152 @@ +#ifndef FiberDRCaloSDAction_H +#define FiberDRCaloSDAction_H + +// DD4hep Framework include files +#include "DD4hep/Version.h" +#include "DD4hep/Segmentations.h" + +#include "DDG4/Geant4Data.h" + +#include + +#include + +#if DD4HEP_VERSION_GE(1, 21) +#define GEANT4_CONST_STEP const +#else +#define GEANT4_CONST_STEP +#endif + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep +{ + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim + { + + class Geant4DRCalorimeter + { + public: + class Hit : public Geant4HitData + { + public: + typedef Geant4HitData base_t; + typedef std::map DRsimTimeStruct; + typedef std::map DRsimWavlenSpectrum; + + /// Hit position + ROOT::Math::XYZVector position; + /// Hit contributions by individual particles + std::vector truth; + /// Total energy deposit + double energyDeposit; + + public: + /// Default constructor (for ROOT) + Hit(); + // Constructor + Hit(float wavSampling, float timeSampling) + : Geant4HitData(), + fSiPMnum(0), + fPhotons(0), + mWavSampling(wavSampling), + mTimeSampling(timeSampling) + { + } + /// copy constructor + Hit(const Hit &right) + : Geant4HitData() + { + fSiPMnum = right.fSiPMnum; + fPhotons = right.fPhotons; + fWavlenSpectrum = right.fWavlenSpectrum; + fTimeStruct = right.fTimeStruct; + mWavSampling = right.mWavSampling; + mTimeSampling = right.mTimeSampling; + } + /// Copy assignment operator + Hit &operator=(Hit &right) + { + fSiPMnum = right.fSiPMnum; + fPhotons = right.fPhotons; + fWavlenSpectrum = right.fWavlenSpectrum; + fTimeStruct = right.fTimeStruct; + mWavSampling = right.mWavSampling; + mTimeSampling = right.mTimeSampling; + return *this; + } + G4bool operator==(const Hit &right) const + { + return (fSiPMnum == right.fSiPMnum); + } + /// Move constructor + Hit(Hit &&c) = delete; + + /// Standard constructor + Hit(const Position &cell_pos); + /// Default destructor + virtual ~Hit(){}; + /// Move assignment operator + Hit &operator=(Hit &&c) = delete; + + void photonCount() { fPhotons++; } + unsigned long GetPhotonCount() const { return fPhotons; } + + void SetSiPMnum(dd4hep::DDSegmentation::CellID n) { fSiPMnum = n; } + const dd4hep::DDSegmentation::CellID &GetSiPMnum() const { return fSiPMnum; } + + void CountWavlenSpectrum(int ibin) + { + auto it = fWavlenSpectrum.find(ibin); + + if (it == fWavlenSpectrum.end()) + fWavlenSpectrum.insert(std::make_pair(ibin, 1)); + else + it->second++; + }; + const DRsimWavlenSpectrum &GetWavlenSpectrum() const { return fWavlenSpectrum; } + + void CountTimeStruct(int ibin) + { + auto it = fTimeStruct.find(ibin); + + if (it == fTimeStruct.end()) + fTimeStruct.insert(std::make_pair(ibin, 1)); + else + it->second++; + }; + const DRsimTimeStruct &GetTimeStruct() const { return fTimeStruct; } + + float GetSamplingTime() const { return mTimeSampling; } + float GetSamplingWavlen() const { return mWavSampling; } + + float GetTimeStart() const { return fTimeStart; } + void SetTimeStart(float val) { fTimeStart = val; } + + float GetTimeEnd() const { return fTimeEnd; } + void SetTimeEnd(float val) { fTimeEnd = val; } + + float GetWavlenMax() const { return fWavlenMax; } + void SetWavlenMax(float val) { fWavlenMax = val; } + + float GetWavlenMin() const { return fWavlenMin; } + void SetWavlenMin(float val) { fWavlenMin = val; } + + private: + dd4hep::DDSegmentation::CellID fSiPMnum; + unsigned long fPhotons; + DRsimWavlenSpectrum fWavlenSpectrum; + DRsimTimeStruct fTimeStruct; + float mWavSampling; + float mTimeSampling; + float fWavlenMax; + float fWavlenMin; + float fTimeStart; + float fTimeEnd; + }; + }; + + } +} +#endif diff --git a/plugins/Geant4Output2EDM4hep_DRC.cpp b/plugins/Geant4Output2EDM4hep_DRC.cpp new file mode 100644 index 000000000..6d408424e --- /dev/null +++ b/plugins/Geant4Output2EDM4hep_DRC.cpp @@ -0,0 +1,695 @@ +#ifndef DD4HEP_DDG4_Geant4Output2EDM4hep_DRC_H +#define DD4HEP_DDG4_Geant4Output2EDM4hep_DRC_H + +/// Framework include files +#include +#include +#include +#include + +/// edm4hep include files +#include +#include +#include +#include +#include +#include +#include +/// podio include files +#include +#include +#include +#if PODIO_BUILD_VERSION >= PODIO_VERSION(0, 99, 0) +#include +#else +#include +namespace podio { + using ROOTWriter = podio::ROOTFrameWriter; +} +#endif + +#include "FiberDRCaloSDAction.h" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + class ComponentCast; + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + class Geant4ParticleMap; + + /// Base class to output Geant4 event data to EDM4hep + /** + * \author F.Gaede + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4Output2EDM4hep_DRC : public Geant4OutputAction { + protected: + using writer_t = podio::ROOTWriter; + using stringmap_t = std::map< std::string, std::string >; + using trackermap_t = std::map< std::string, edm4hep::SimTrackerHitCollection >; + using calorimeterpair_t = std::pair< edm4hep::SimCalorimeterHitCollection, edm4hep::CaloHitContributionCollection >; + using calorimetermap_t = std::map< std::string, calorimeterpair_t >; + using drcalopair_t = std::pair< edm4hep::RawCalorimeterHitCollection, edm4hep::RawTimeSeriesCollection> ; // Required info for IDEA DRC sim hit + using drcalomap_t = std::map< std::string, drcalopair_t >; // Required info for IDEA DRC sim hit + using drcaloWavmap_t = std::map< std::string, edm4hep::RawTimeSeriesCollection >; + + std::unique_ptr m_file { }; + podio::Frame m_frame { }; + edm4hep::MCParticleCollection m_particles { }; + trackermap_t m_trackerHits; + calorimetermap_t m_calorimeterHits; + drcalomap_t m_drcaloHits; + drcaloWavmap_t m_drcaloWaves; + stringmap_t m_runHeader; + stringmap_t m_eventParametersInt; + stringmap_t m_eventParametersFloat; + stringmap_t m_eventParametersString; + stringmap_t m_cellIDEncodingStrings{}; + std::string m_section_name { "events" }; + int m_runNo { 0 }; + int m_runNumberOffset { 0 }; + int m_eventNo { 0 }; + int m_eventNumberOffset { 0 }; + bool m_filesByRun { false }; + + /// Data conversion interface for MC particles to EDM4hep format + void saveParticles(Geant4ParticleMap* particles); + /// Store the metadata frame with e.g. the cellID encoding strings + void saveFileMetaData(); + public: + /// Standard constructor + Geant4Output2EDM4hep_DRC(Geant4Context* ctxt, const std::string& nam); + /// Default destructor + virtual ~Geant4Output2EDM4hep_DRC(); + /// Callback to store the Geant4 run information + virtual void beginRun(const G4Run* run); + /// Callback to store the Geant4 run information + virtual void endRun(const G4Run* run); + + /// Callback to store the Geant4 run information + virtual void saveRun(const G4Run* run); + /// Callback to store the Geant4 event + virtual void saveEvent( OutputContext& ctxt); + /// Callback to store each Geant4 hit collection + virtual void saveCollection( OutputContext& ctxt, G4VHitsCollection* collection); + /// Commit data at end of filling procedure + virtual void commit( OutputContext& ctxt); + + /// begin-of-event callback - creates EDM4hep event and adds it to the event context + virtual void begin(const G4Event* event); + protected: + /// Fill event parameters in EDM4hep event + template + void saveEventParameters(const std::map& parameters) { + for(const auto& p : parameters) { + info("Saving event parameter: %-32s = %s", p.first.c_str(), p.second.c_str()); + m_frame.putParameter(p.first, p.second); + } + } + }; + + template <> void EventParameters::extractParameters(podio::Frame& frame) { + for(auto const& p: this->intParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } + for(auto const& p: this->fltParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } + for(auto const& p: this->strParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } +#if PODIO_BUILD_VERSION > PODIO_VERSION(0, 16, 2) + // This functionality is only present in podio > 0.16.2 + for (auto const& p: this->dblParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } +#endif + } + + template <> void RunParameters::extractParameters(podio::Frame& frame) { + for(auto const& p: this->intParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } + for(auto const& p: this->fltParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } + for(auto const& p: this->strParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } +#if PODIO_BUILD_VERSION > PODIO_VERSION(0, 16, 2) + // This functionality is only present in podio > 0.16.2 + for (auto const& p: this->dblParameters()) { + printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str()); + frame.putParameter(p.first, p.second); + } +#endif + } + + } // End namespace sim +} // End namespace dd4hep +#endif // DD4HEP_DDG4_Geant4Output2EDM4hep_DRC_H + +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : F.Gaede, DESY +// +//========================================================================== + +/// Framework include files +#include +#include + +#include +#include +#include +#include +#include +#include + +///#include +/// Geant4 headers +#include +#include +#include +#include +#include +#include +#include +/// use the Geant4 units in namespace CLHEP +#include + +/// edm4hep include files +#include + +using namespace dd4hep::sim; +using namespace dd4hep; + +namespace { + G4Mutex action_mutex = G4MUTEX_INITIALIZER; +} + +#include +DECLARE_GEANT4ACTION(Geant4Output2EDM4hep_DRC) + +/// Standard constructor +Geant4Output2EDM4hep_DRC::Geant4Output2EDM4hep_DRC(Geant4Context* ctxt, const std::string& nam) +: Geant4OutputAction(ctxt,nam), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0) +{ + declareProperty("RunHeader", m_runHeader); + declareProperty("EventParametersInt", m_eventParametersInt); + declareProperty("EventParametersFloat", m_eventParametersFloat); + declareProperty("EventParametersString", m_eventParametersString); + declareProperty("RunNumberOffset", m_runNumberOffset); + declareProperty("EventNumberOffset", m_eventNumberOffset); + declareProperty("SectionName", m_section_name); + declareProperty("FilesByRun", m_filesByRun); + info("Writer is now instantiated ..." ); + InstanceCount::increment(this); +} + +/// Default destructor +Geant4Output2EDM4hep_DRC::~Geant4Output2EDM4hep_DRC() { + G4AutoLock protection_lock(&action_mutex); + m_file.reset(); + InstanceCount::decrement(this); +} + +// Callback to store the Geant4 run information +void Geant4Output2EDM4hep_DRC::beginRun(const G4Run* run) { + G4AutoLock protection_lock(&action_mutex); + std::string fname = m_output; + m_runNo = run->GetRunID(); + if ( m_filesByRun ) { + std::size_t idx = m_output.rfind("."); + if ( idx != std::string::npos ) { + fname = m_output.substr(0, idx) + _toString(m_runNo, ".run%08d") + m_output.substr(idx); + } + } + if ( !fname.empty() ) { + m_file = std::make_unique(fname); + if ( !m_file ) { + fatal("+++ Failed to open output file: %s", fname.c_str()); + } + printout( INFO, "Geant4Output2EDM4hep_DRC" ,"Opened %s for output", fname.c_str() ) ; + } +} + +/// Callback to store the Geant4 run information +void Geant4Output2EDM4hep_DRC::endRun(const G4Run* run) { + saveRun(run); + saveFileMetaData(); + if ( m_file ) { + m_file->finish(); + m_file.reset(); + } +} + +void Geant4Output2EDM4hep_DRC::saveFileMetaData() { + podio::Frame metaFrame{}; + for (const auto& [name, encodingStr] : m_cellIDEncodingStrings) { + metaFrame.putParameter(name + "__CellIDEncoding", encodingStr); + } + + m_file->writeFrame(metaFrame, "metadata"); +} + +/// Commit data at end of filling procedure +void Geant4Output2EDM4hep_DRC::commit( OutputContext& /* ctxt */) { + if ( m_file ) { + G4AutoLock protection_lock(&action_mutex); + m_frame.put( std::move(m_particles), "MCParticles"); + for (auto it = m_trackerHits.begin(); it != m_trackerHits.end(); ++it) { + m_frame.put( std::move(it->second), it->first); + } + for (auto& [colName, calorimeterHits] : m_calorimeterHits) { + m_frame.put( std::move(calorimeterHits.first), colName); + m_frame.put( std::move(calorimeterHits.second), colName + "Contributions"); + } + for (auto& [colName, calorimeterHits] : m_drcaloHits) { + m_frame.put( std::move(calorimeterHits.first), colName + "RawHit"); + m_frame.put( std::move(calorimeterHits.second), colName + "TimeStruct"); + } + for (auto it = m_drcaloWaves.begin(); it != m_drcaloWaves.end(); ++it) { + m_frame.put( std::move(it->second), it->first + "WaveLen"); + } + m_file->writeFrame(m_frame, m_section_name); + m_particles.clear(); + m_trackerHits.clear(); + m_calorimeterHits.clear(); + m_drcaloHits.clear(); + m_drcaloWaves.clear(); + m_frame = {}; + return; + } + except("+++ Failed to write output file. [Stream is not open]"); +} + +/// Callback to store the Geant4 run information +void Geant4Output2EDM4hep_DRC::saveRun(const G4Run* run) { + G4AutoLock protection_lock(&action_mutex); + // --- write an edm4hep::RunHeader --------- + // Runs are just Frames with different contents in EDM4hep / podio. We simply + // store everything as parameters for now + podio::Frame runHeader {}; + for (const auto& [key, value] : m_runHeader) + runHeader.putParameter(key, value); + + m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID(); + runHeader.putParameter("runNumber", m_runNo); + runHeader.putParameter("GEANT4Version", G4Version); + runHeader.putParameter("DD4hepVersion", versionString()); + runHeader.putParameter("detectorName", context()->detectorDescription().header().name()); + + RunParameters* parameters = context()->run().extension(false); + if ( parameters ) { + parameters->extractParameters(runHeader); + } + + m_file->writeFrame(runHeader, "runs"); +} + +void Geant4Output2EDM4hep_DRC::begin(const G4Event* event) { + /// Create event frame object + m_eventNo = event->GetEventID(); + m_frame = {}; + m_particles = {}; + m_trackerHits.clear(); + m_calorimeterHits.clear(); + m_drcaloHits.clear(); + m_drcaloWaves.clear(); +} + +/// Data conversion interface for MC particles to EDM4hep format +void Geant4Output2EDM4hep_DRC::saveParticles(Geant4ParticleMap* particles) { + typedef detail::ReferenceBitMask PropertyMask; + typedef Geant4ParticleMap::ParticleMap ParticleMap; + const ParticleMap& pm = particles->particleMap; + + m_particles.clear(); + if ( pm.size() > 0 ) { + size_t cnt = 0; + // Mapping of ids in the ParticleMap to indices in the MCParticle collection + std::map p_ids; + std::vector p_part; + p_part.reserve(pm.size()); + // First create the particles + for (const auto& iParticle : pm) { + int id = iParticle.first; + const Geant4ParticleHandle p = iParticle.second; + PropertyMask mask(p->status); + // std::cout << " ********** mcp status : 0x" << std::hex << p->status << ", mask.isSet(G4PARTICLE_GEN_STABLE) x" << std::dec << mask.isSet(G4PARTICLE_GEN_STABLE) <pdgID); + // Because EDM4hep is switching between vector3f[loat] and vector3d[ouble] + using MT = decltype(std::declval().getMomentum().x); + mcp.setMomentum( {MT(p->psx/CLHEP::GeV),MT(p->psy/CLHEP::GeV),MT(p->psz/CLHEP::GeV)} ); + mcp.setMomentumAtEndpoint( {MT(p->pex/CLHEP::GeV),MT(p->pey/CLHEP::GeV),MT(p->pez/CLHEP::GeV)} ); + + double vs_fa[3] = { p->vsx/CLHEP::mm, p->vsy/CLHEP::mm, p->vsz/CLHEP::mm } ; + mcp.setVertex( vs_fa ); + + double ve_fa[3] = { p->vex/CLHEP::mm, p->vey/CLHEP::mm, p->vez/CLHEP::mm } ; + mcp.setEndpoint( ve_fa ); + + mcp.setTime(p->time/CLHEP::ns); + mcp.setMass(p->mass/CLHEP::GeV); + mcp.setCharge(def ? def->GetPDGCharge() : 0); // Charge(e+) = 1 ! + + // Set generator status + mcp.setGeneratorStatus(0); + if( p->genStatus ) { + mcp.setGeneratorStatus( p->genStatus ) ; + } else { + if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) mcp.setGeneratorStatus(1); + else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) mcp.setGeneratorStatus(2); + else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3); + else if ( mask.isSet(G4PARTICLE_GEN_BEAM) ) mcp.setGeneratorStatus(4); + else if ( mask.isSet(G4PARTICLE_GEN_OTHER) ) mcp.setGeneratorStatus(9); + } + + // Set simulation status + mcp.setCreatedInSimulation( mask.isSet(G4PARTICLE_SIM_CREATED) ); + mcp.setBackscatter( mask.isSet(G4PARTICLE_SIM_BACKSCATTER) ); + mcp.setVertexIsNotEndpointOfParent( mask.isSet(G4PARTICLE_SIM_PARENT_RADIATED) ); + mcp.setDecayedInTracker( mask.isSet(G4PARTICLE_SIM_DECAY_TRACKER) ); + mcp.setDecayedInCalorimeter( mask.isSet(G4PARTICLE_SIM_DECAY_CALO) ); + mcp.setHasLeftDetector( mask.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) ); + mcp.setStopped( mask.isSet(G4PARTICLE_SIM_STOPPED) ); + mcp.setOverlay( false ); + + //fg: if simstatus !=0 we have to set the generator status to 0: + if( mcp.isCreatedInSimulation() ) + mcp.setGeneratorStatus( 0 ) ; + + mcp.setSpin(p->spin); + mcp.setColorFlow(p->colorFlow); + + p_ids[id] = cnt++; + p_part.push_back(p); + } + + // Now establish parent-daughter relationships + for(size_t i=0; i < p_ids.size(); ++i) { + const Geant4Particle* p = p_part[i]; + auto q = m_particles[i]; + + for (const auto& idau : p->daughters) { + const auto k = p_ids.find(idau); + if (k == p_ids.end()) { + fatal("+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau); + continue; + } + int iqdau = (*k).second; + auto qdau = m_particles[iqdau]; + q.addToDaughters(qdau); + } + + for (const auto& ipar : p->parents) { + if (ipar >= 0) { // A parent ID of -1 means NO parent, because a base of 0 is perfectly legal + const auto k = p_ids.find(ipar); + if (k == p_ids.end()) { + fatal("+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar); + continue; + } + int iqpar = (*k).second; + auto qpar = m_particles[iqpar]; + q.addToParents(qpar); + } + } + } + } +} + +/// Callback to store the Geant4 event +void Geant4Output2EDM4hep_DRC::saveEvent(OutputContext& ctxt) { + EventParameters* parameters = context()->event().extension(false); + int runNumber(0), eventNumber(0); + const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0); + const int runNumberOffset(m_runNumberOffset > 0 ? m_runNumberOffset : 0); + double eventWeight{0}; + // Get event number, run number and parameters from extension ... + if ( parameters ) { + runNumber = parameters->runNumber() + runNumberOffset; + eventNumber = parameters->eventNumber() + eventNumberOffset; + parameters->extractParameters(m_frame); +#if PODIO_BUILD_VERSION > PODIO_VERSION(0, 16, 2) + // This functionality is only present in podio > 0.16.2 +#if PODIO_BUILD_VERSION > PODIO_VERSION(0, 99, 0) + eventWeight = m_frame.getParameter("EventWeights").value_or(0.0); +#else + eventWeight = m_frame.getParameter("EventWeights"); +#endif +#endif + } else { // ... or from DD4hep framework + runNumber = m_runNo + runNumberOffset; + eventNumber = ctxt.context->GetEventID() + eventNumberOffset; + } + printout(INFO,"Geant4Output2EDM4hep_DRC","+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber); + + // this does not compile as create() is we only get a const ref - need to review PODIO EventStore API + edm4hep::EventHeaderCollection header_collection; + auto header = header_collection.create(); + header.setRunNumber(runNumber); + header.setEventNumber(eventNumber); + header.setWeight(eventWeight); + //not implemented in EDM4hep ? header.setDetectorName(context()->detectorDescription().header().name()); + header.setTimeStamp( std::time(nullptr) ) ; + m_frame.put( std::move(header_collection), "EventHeader"); + + saveEventParameters(m_eventParametersInt); + saveEventParameters(m_eventParametersFloat); + saveEventParameters(m_eventParametersString); + + Geant4ParticleMap* part_map = context()->event().extension(false); + if ( part_map ) { + print("+++ Saving %d EDM4hep particles....",int(part_map->particleMap.size())); + if ( part_map->particleMap.size() > 0 ) { + saveParticles(part_map); + } + } +} + +/** + * Helper struct that can be used together with map::try_emplace to construct + * the encoding only once per collection (name). + */ +struct LazyEncodingExtraction { + /// Constructor that does effectively nothing. This will be called in every + /// try_emplace call + LazyEncodingExtraction(Geant4HitCollection* coll) : m_coll(coll) {} + /// Defer the real work to the implicit conversion to std::string that will + /// only be called if the value is actually emplaced into the map + operator std::string() const { + const auto* sd = m_coll->sensitive(); + return dd4hep::sim::Geant4ConversionHelper::encoding(sd->sensitiveDetector()); + } +private: + Geant4HitCollection* m_coll{nullptr}; +}; + + +/// Callback to store each Geant4 hit collection +void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/, G4VHitsCollection* collection) { + Geant4HitCollection* coll = dynamic_cast(collection); + std::string colName = collection->GetName(); + if( coll == nullptr ){ + error(" no Geant4HitCollection: %s ", colName.c_str()); + return ; + } + size_t nhits = collection->GetSize(); + Geant4ParticleMap* pm = context()->event().extension(false); + debug("+++ Saving EDM4hep collection %s with %d entries.", colName.c_str(), int(nhits)); + + // Using try_emplace here to only fill this the first time we come across + m_cellIDEncodingStrings.try_emplace(colName, LazyEncodingExtraction{coll}); + + //------------------------------------------------------------------- + if( typeid( Geant4Tracker::Hit ) == coll->type().type() ){ + // Create the hit container even if there are no entries! + auto& hits = m_trackerHits[colName]; + for(unsigned i=0 ; i < nhits ; ++i){ + auto sth = hits->create(); + const Geant4Tracker::Hit* hit = coll->hit(i); + const Geant4Tracker::Hit::Contribution& t = hit->truth; + int trackID = pm->particleID(t.trackID); + auto mcp = m_particles.at(trackID); + const auto& mom = hit->momentum; + const auto& pos = hit->position; + edm4hep::Vector3f(); + sth.setCellID( hit->cellID ) ; + sth.setEDep(hit->energyDeposit/CLHEP::GeV); + sth.setPathLength(hit->length/CLHEP::mm); + sth.setTime(hit->truth.time/CLHEP::ns); +#if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 10, 99) + sth.setParticle(mcp); +#else + sth.setMCParticle(mcp); +#endif + sth.setPosition( {pos.x()/CLHEP::mm, pos.y()/CLHEP::mm, pos.z()/CLHEP::mm} ); + sth.setMomentum( {float(mom.x()/CLHEP::GeV),float(mom.y()/CLHEP::GeV),float(mom.z()/CLHEP::GeV)} ); + auto particleIt = pm->particles().find(trackID); + if( ( particleIt != pm->particles().end()) ){ + // if the original track ID of the particle is not the same as the + // original track ID of the hit it was produced by an MCParticle that + // is no longer stored + sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.trackID) ); + } + } + //------------------------------------------------------------------- + } + else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type() ){ + Geant4Sensitive* sd = coll->sensitive(); + int hit_creation_mode = sd->hitCreationMode(); + // Create the hit container even if there are no entries! + auto& hits = m_calorimeterHits[colName]; + for(unsigned i=0 ; i < nhits ; ++i){ + auto sch = hits.first->create(); + const Geant4Calorimeter::Hit* hit = coll->hit(i); + const auto& pos = hit->position; + sch.setCellID( hit->cellID ); + sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)}); + sch.setEnergy( hit->energyDeposit/CLHEP::GeV ); + + + // now add the individual step contributions + for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){ + + auto sCaloHitCont = hits.second->create(); + sch.addToContributions( sCaloHitCont ); + + const Geant4HitData::Contribution& c = *ci; + int trackID = pm->particleID(c.trackID); + auto mcp = m_particles.at(trackID); + sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV ); + sCaloHitCont.setTime( c.time/CLHEP::ns ); + sCaloHitCont.setParticle( mcp ); + + if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE ) { + edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm); + sCaloHitCont.setPDG( c.pdgID ); + sCaloHitCont.setStepPosition( p ); + } + } + } + //------------------------------------------------------------------- + } + else if( typeid( Geant4DRCalorimeter::Hit ) == coll->type().type() ){ + Geant4Sensitive* sd = coll->sensitive(); + int hit_creation_mode = sd->hitCreationMode(); + // Create the hit container even if there are no entries! + auto& hits = m_calorimeterHits[colName]; + auto& DRhits = m_drcaloHits[colName]; + auto& DRwaves = m_drcaloWaves[colName]; + + for(unsigned i=0 ; i < nhits ; ++i){ + auto sch = hits.first->create(); + const Geant4DRCalorimeter::Hit* hit = coll->hit(i); + const auto& pos = hit->position; + sch.setCellID( hit->cellID ); + sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)}); + sch.setEnergy( hit->energyDeposit/CLHEP::GeV ); + + // now add the individual step contributions + for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){ + + auto sCaloHitCont = hits.second->create(); + sch.addToContributions( sCaloHitCont ); + + const Geant4HitData::Contribution& c = *ci; + int trackID = pm->particleID(c.trackID); + auto mcp = m_particles.at(trackID); + sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV ); + sCaloHitCont.setTime( c.time/CLHEP::ns ); + sCaloHitCont.setParticle( mcp ); + + if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE ) { + edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm); + sCaloHitCont.setPDG( c.pdgID ); + sCaloHitCont.setStepPosition( p ); + } + } + + // For DRC raw calo hit & raw time series + auto rawCaloHits = DRhits.first->create(); + auto rawTimeStruct = DRhits.second->create(); + auto rawWaveStruct = DRwaves->create(); + + float samplingT = hit->GetSamplingTime(); + float timeStart = hit->GetTimeStart(); + float timeEnd = hit->GetTimeEnd(); + auto& timemap = hit->GetTimeStruct(); + + rawTimeStruct.setInterval(samplingT); + rawTimeStruct.setTime(timeStart); + rawTimeStruct.setCharge( static_cast(hit->GetPhotonCount()) ); + rawTimeStruct.setCellID( hit->cellID ); + + // abuse time series for saving wavelength spectrum (for R&D purpose) + float samplingW = hit->GetSamplingWavlen(); + float wavMax = hit->GetWavlenMax(); + float wavMin = hit->GetWavlenMin(); + auto& wavmap = hit->GetWavlenSpectrum(); + rawWaveStruct.setInterval(samplingW); + rawWaveStruct.setTime(wavMin); + rawWaveStruct.setCharge( static_cast(hit->GetPhotonCount()) ); + rawWaveStruct.setCellID( hit->cellID ); + + unsigned nbinTime = static_cast((timeEnd-timeStart)/samplingT); + unsigned nbinWav = static_cast((wavMax-wavMin)/samplingW); + int peakTime = 0.; + int peakVal = 0; + + // same as the ROOT TH1 binning scheme (0: underflow, nbin+1:overflow) + for (unsigned itime = 1; itime < nbinTime+1; itime++) { + int count = 0; + + if ( timemap.find(itime)!=timemap.end() ) + count = timemap.at(itime); + + int candidate = std::max( peakVal, count ); + + if ( peakVal < candidate ) { + peakVal = candidate; + peakTime = itime; + } + + rawTimeStruct.addToAdcCounts(count); + } + + for (unsigned iwav = 1; iwav < nbinWav+1; iwav++) { + int count = 0; + + if ( wavmap.find(iwav)!=wavmap.end() ) + count = wavmap.at(iwav); + + rawWaveStruct.addToAdcCounts(count); + } + + rawCaloHits.setCellID( hit->cellID ); + rawCaloHits.setAmplitude( hit->GetPhotonCount() ); + rawCaloHits.setTimeStamp( peakTime-1 + static_cast(timeStart/samplingT) ); + } + //------------------------------------------------------------------- + } else { + error("+++ unknown type in Geant4HitCollection %s ", coll->type().type().name()); + } +} \ No newline at end of file diff --git a/plugins/LinearSortingPolicy.cpp b/plugins/LinearSortingPolicy.cpp index dde46cf9a..a84dda5e0 100644 --- a/plugins/LinearSortingPolicy.cpp +++ b/plugins/LinearSortingPolicy.cpp @@ -26,7 +26,6 @@ #include #include -#include #include using dd4hep::DetElement; diff --git a/plugins/TPCSDAction.cpp b/plugins/TPCSDAction.cpp index 5123eff09..b9f4322db 100644 --- a/plugins/TPCSDAction.cpp +++ b/plugins/TPCSDAction.cpp @@ -2,7 +2,6 @@ #include "DDG4/Geant4SensDetAction.inl" #include "DDG4/Geant4EventAction.h" #include "DDG4/Geant4Mapping.h" -#include "G4OpticalPhoton.hh" #include "G4VProcess.hh" #if DD4HEP_VERSION_GE(1, 21) @@ -453,4 +452,5 @@ namespace dd4hep { #include "DDG4/Factories.h" + DECLARE_GEANT4SENSITIVE( TPCSDAction ) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3d4876777..5ba8150b8 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -91,10 +91,20 @@ SET_TESTS_PROPERTIES( t_${test_name} PROPERTIES FAIL_REGULAR_EXPRESSION "Except if(DCH_INFO_H_EXIST) SET( test_name "test_IDEA_o1_v03" ) ADD_TEST( t_${test_name} "${CMAKE_INSTALL_PREFIX}/bin/run_test_${PackageName}.sh" - ddsim --compactFile=${CMAKE_CURRENT_SOURCE_DIR}/../FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml --runType=batch -G -N=1 --outputFile=testIDEA_o1_v03.slcio ) -SET_TESTS_PROPERTIES( t_${test_name} PROPERTIES FAIL_REGULAR_EXPRESSION " Exception; EXCEPTION;ERROR;Error" TIMEOUT 700) + ddsim --compactFile=${CMAKE_CURRENT_SOURCE_DIR}/../FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml -N 1 -G --gun.distribution uniform --random.seed 1988301045 ) + SET_TESTS_PROPERTIES( t_${test_name} PROPERTIES FAIL_REGULAR_EXPRESSION " Exception; EXCEPTION;ERROR;Error" TIMEOUT 100) endif() +#-------------------------------------------------- +# test for IDEA o1 v03, with DRC +if(DCH_INFO_H_EXIST) +SET( test_name "test_IDEA_with_DRC_o1_v03" ) +ADD_TEST( t_${test_name} "${CMAKE_INSTALL_PREFIX}/bin/run_test_${PackageName}.sh" + ddsim --compactFile=${CMAKE_CURRENT_SOURCE_DIR}/compact/IDEA_withDRC_o1_v03.xml --steeringFile=${CMAKE_CURRENT_SOURCE_DIR}/../example/SteeringFile_IDEA_o1_v03.py -G --gun.distribution uniform --random.seed 1988301045 ) + SET_TESTS_PROPERTIES( t_${test_name} PROPERTIES FAIL_REGULAR_EXPRESSION " Exception; EXCEPTION;ERROR;Error" TIMEOUT 600) +endif() + + #-------------------------------------------------- # test for ALLEGRO o1 v02 SET( test_name "test_ALLEGRO_o1_v02" ) @@ -149,3 +159,13 @@ ADD_TEST( t_SensThickness_Clic_o2_v4 "${CMAKE_INSTALL_PREFIX}/bin/run_test_${Pac ${CMAKE_INSTALL_PREFIX}/bin/TestSensThickness ${CMAKE_CURRENT_SOURCE_DIR}/../CLIC/compact/CLIC_o2_v04/CLIC_o2_v04.xml 300 50 ) ADD_TEST( t_SensThickness_CLIC_o3_v15 "${CMAKE_INSTALL_PREFIX}/bin/run_test_${PackageName}.sh" ${CMAKE_INSTALL_PREFIX}/bin/TestSensThickness ${CMAKE_CURRENT_SOURCE_DIR}/../CLIC/compact/CLIC_o3_v15/CLIC_o3_v15.xml 100 50 ) + +#-------------------------------------------------- +# check if files named the same contain the same in FCCee +ADD_TEST( + NAME t_test_files_versions_FCCee + COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_${PackageName}.sh" + "${PROJECT_SOURCE_DIR}/utils/compareIdenticalFiles.sh" + "${PROJECT_SOURCE_DIR}/FCCee" + "${PROJECT_SOURCE_DIR}/utils/IdenticalFiles_ToBeIgnored.txt" +) diff --git a/test/compact/IDEA_withDRC_o1_v03.xml b/test/compact/IDEA_withDRC_o1_v03.xml new file mode 100644 index 000000000..22ab501bb --- /dev/null +++ b/test/compact/IDEA_withDRC_o1_v03.xml @@ -0,0 +1,99 @@ + + + + + + + + Version o1_v03 of the IDEA detector + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/src/BeamCalZtest.cpp b/test/src/BeamCalZtest.cpp index 4d6760806..ea08ac0cc 100644 --- a/test/src/BeamCalZtest.cpp +++ b/test/src/BeamCalZtest.cpp @@ -1,6 +1,5 @@ #include #include -#include "DD4hep/DD4hepUnits.h" #include #include diff --git a/test/src/TestSensThickness.cpp b/test/src/TestSensThickness.cpp index a777a3631..9f6162fd2 100644 --- a/test/src/TestSensThickness.cpp +++ b/test/src/TestSensThickness.cpp @@ -1,14 +1,9 @@ // Test the setting of the sensitiveThickness for the DataStructs in some Tracker detectorDrivers -#include #include #include -#include #include -#include -#include -#include #include #include #include diff --git a/utils/IdenticalFiles_ToBeIgnored.txt b/utils/IdenticalFiles_ToBeIgnored.txt new file mode 100644 index 000000000..9e760a722 --- /dev/null +++ b/utils/IdenticalFiles_ToBeIgnored.txt @@ -0,0 +1,35 @@ +README.md +CMakeLists.txt +DectDimensions.xml +materials.xml +#__CLD_ill_named_files__ +additional_materials.xml +BeamInstrumentation_o3_v01_overlap.xml +BeamInstrumentation_o3_v02_fitShield.xml +Beampipe_o4_v03_noNotch_Ta_cone.xml +Beampipe_o4_v04_noNotch_W_n02_smallBP.xml +Beampipe_o4_v04_noNotch_W_n02.xml +Beampipe_o4_v05.xml +ECalBarrel_o2_v01_03.xml +ECalEndcap_o2_v01_03.xml +elements.xml +FCCee_DectDimensions.xml +FCCee_o2_v02.xml +HCalBarrel_o1_v01_01.xml +HCalEndcap_o1_v01_01.xml +HOMAbsorber_v00.xml +InnerTracker_o2_v06_02.xml +LumiCal_o3_v02_02.xml +LumiCal_o3_v02_03.xml +OuterTracker_o2_v06_02.xml +Solenoid_o1_v01_02.xml +Vertex_o4_v05_smallBP.xml +Vertex_o4_v05.xml +YokeBarrel_o1_v01_02.xml +YokeEndcap_o1_v01_02.xml +#__ALLEGRO_ill_named_files__ +ECalBarrel_thetamodulemerged.xml +ECalBarrel_thetamodulemerged_calibration.xml +ECalBarrel_thetamodulemerged_upstream.xml +#__FCCee_ill_named_files__ +FCCee_DectEmptyMaster.xml diff --git a/utils/compareIdenticalFiles.sh b/utils/compareIdenticalFiles.sh new file mode 100755 index 000000000..e22c658c0 --- /dev/null +++ b/utils/compareIdenticalFiles.sh @@ -0,0 +1,49 @@ +#!/bin/bash + +# This script compares all the files named the same within a given path (as argument) +# If files with the same name have different contents, it prints the paths of those files +# If the files are identical, nothing is printed +# Files listed in an ignore file are excluded + +# Usage: ./script.sh + +search_dir="$1" +ignore_file="$2" + +# Read filenames to ignore from the ignore file +files_to_ignore=$(<"$ignore_file") + +# Create an associative array to store file paths with the same names +declare -A file_names + +# Initialize status code variable +status_code=0 + +# Iterate through files in the directory +while IFS= read -r -d '' file; do + # Get the base name of the file (without the path) + file_name=$(basename "$file") + + # Ignore files listed in the ignore file + if grep -qx "$file_name" <<< "$files_to_ignore"; then + continue + fi + + # Check if there is already a file with the same name + if [[ -n "${file_names[$file_name]}" ]]; then + # Compare the contents of the files + if ! cmp -s "$file" "${file_names[$file_name]}"; then + # Print the paths of the files if they differ + echo "Error. Files with the same name but different contents:" + echo "$file" + echo "${file_names[$file_name]}" + echo + status_code=1 + fi + else + # Store the file path in the array + file_names["$file_name"]="$file" + fi +done < <(find "$search_dir" -not -path '*/.*' -type f -print0) + +exit $status_code