diff --git a/.github/scripts/materialScanAxis.sh b/.github/scripts/materialScanAxis.sh new file mode 100755 index 00000000..c3d2661c --- /dev/null +++ b/.github/scripts/materialScanAxis.sh @@ -0,0 +1,57 @@ +#!/bin/bash + +# This test scan the material along the X (or Z) axis using materialScan from DD4hep +# then extract the second-to-last line, which contains a summary of the scan, +# then checks if X0 is larger or smaller than X0 for air + +compactFile=$1 +axis=$2 +IsMaterialExpectedToBeAir=$( awk '{print tolower($0)}' <<< $3) + + +#Check if file exists... +[ ! -f "$compactFile" ] && { echo "Compact file $compactFile not found."; exit 2; } + + +#Check if materialScan command from DD4hep if found +if ! command -v materialScan &> /dev/null +then + echo "materialScan could not be found" + exit 3 +fi + + + +#Get the total X0 integrated along the axis. X0 is provided as 12th parameter of the second-to-last line of +x0="" +if [ "x" = $axis ] || [ "X" = $axis ] ; then + x0=$(materialScan $compactFile 0 0 0 1000 0 0 | awk -F' ' '{prevlast = last; last = $0} END {if (NR >= 2) print prevlast}' | awk -F' ' '{ print $12 }') + +elif [ "z" = $axis ] || [ "Z" = $axis ] ; then + x0=$(materialScan $compactFile 0 0 0 0 0 1000 | awk -F' ' '{prevlast = last; last = $0} END {if (NR >= 2) print prevlast}' | head -n 1 | awk -F' ' '{ print $12 }') + +else + echo "Axis " $axis " not supported" + exit 22 +fi + +x0=$(expr "$x0" ) +echo "x0 is " $x0 "cm" + + +xair=0.032756 +if [ "air" = $IsMaterialExpectedToBeAir ] ; then + if (( $(echo "$x0 > $xair" |bc -l) )); then + echo "X0 larger than or equal to air X0 (0.032756 cm)" + exit 22 + else + exit 0 + fi +else + if (( $(echo "$x0 <= $xair" |bc -l) )); then + echo "X0 smaller than or equal to air X0 (0.032756 cm)" + exit 22 + else + exit 0 + fi +fi \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index ffd123d0..ba27ec6a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,7 +31,7 @@ find_package(DD4hep) get_target_property(ddcore_lib DD4hep::DDCore LOCATION) get_filename_component(ddcore_loc ${ddcore_lib} DIRECTORY) function(k4_set_test_env _testname) -foreach(_subdir DetCommon DetFCCeeCLD DetFCCeeCommon DetFCCeeECalInclined DetFCCeeHCalTile DetFCCeeIDEA DetFCCeeIDEA-LAr DetFCChhBaseline1 DetFCChhCalDiscs DetFCChhECalInclined DetFCChhECalSimple DetFCChhHCalTile DetFCChhTailCatcher DetFCChhTrackerTkLayout DetSegmentation DetSensitive DetTrackerSimple) +foreach(_subdir DetCommon DetFCCeeCLD DetFCCeeCommon DetFCCeeECalInclined DetFCCeeHCalTile DetFCCeeIDEA DetFCCeeIDEA-LAr DetFCChhBaseline1 DetFCChhCalDiscs DetFCChhECalInclined DetFCChhECalSimple DetFCChhHCalTile DetFCChhTailCatcher DetFCChhTrackerTkLayout DetSegmentation DetSensitive DetTrackerSimple DRcalo DRcomponents DRsegmentation DRsensitive ) set( _subdirs "${CMAKE_BINARY_DIR}/Detector/${_subdir}:${_subdirs}") endforeach() set_property(TEST ${_testname} APPEND PROPERTY ENVIRONMENT "LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:${CMAKE_CURRENT_BINARY_DIR}:${ddcore_loc}:${_subdirs}:$ENV{LD_LIBRARY_PATH}") diff --git a/Detector/CMakeLists.txt b/Detector/CMakeLists.txt index 7679982f..0b94318f 100644 --- a/Detector/CMakeLists.txt +++ b/Detector/CMakeLists.txt @@ -18,4 +18,12 @@ add_subdirectory(DetFCCeeCLD) add_subdirectory(DetFCCeeCommon) add_subdirectory(DetFCCeeHCalTile) add_subdirectory(DetFCCeeIDEA) + + +add_subdirectory(DRsegmentation) +add_subdirectory(DRcomponents) +add_subdirectory(DRcalo) +add_subdirectory(DRsensitive) + add_subdirectory(DetFCCeeCalDiscs) + diff --git a/Detector/DRcalo/CMakeLists.txt b/Detector/DRcalo/CMakeLists.txt new file mode 100644 index 00000000..baf5d9f7 --- /dev/null +++ b/Detector/DRcalo/CMakeLists.txt @@ -0,0 +1,46 @@ +################################################################################ +# Package: Dual Readout Calorimeter +################################################################################ + +file(GLOB headers ${CMAKE_CURRENT_SOURCE_DIR}/include/*.h +) + +dd4hep_add_plugin(ddDRcalo SOURCES src/*.cpp USES + DD4hep::DDCore + DD4hep::DDCond + ROOT::Core + ROOT::Geom + ROOT::GenVector + ROOT::MathCore + DRsegmentation +) + +target_include_directories(ddDRcalo PUBLIC + $ + $ +) + +set_target_properties(ddDRcalo PROPERTIES PUBLIC_HEADER "${headers}") + +install(TARGETS ddDRcalo + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" COMPONENT dev +) + +install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/compact DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}/Detector/DRcalo) + +# This test scan the material along the X (or Z) axis, +# then extract the second-to-last line, which contains a summary of the scan, +# then fail if the argument 12, which correspond to the integrated X0 is less than the X0 of air + +add_test(NAME ScanMaterialDRcaloX + COMMAND bash ".github/scripts/materialScanAxis.sh" "Detector/DRcalo/compact/DRcalo.xml" "x" "noair" + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} + ) +add_test(NAME ScanMaterialDRcaloZ + COMMAND bash ".github/scripts/materialScanAxis.sh" "Detector/DRcalo/compact/DRcalo.xml" "z" "air" + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} + ) + +k4_set_test_env(ScanMaterialDRcaloX) +k4_set_test_env(ScanMaterialDRcaloZ) diff --git a/Detector/DRcalo/compact/DRcalo.xml b/Detector/DRcalo/compact/DRcalo.xml new file mode 100644 index 00000000..716fe005 --- /dev/null +++ b/Detector/DRcalo/compact/DRcalo.xml @@ -0,0 +1,745 @@ + + + + + + + + + + + + The compact format for the dual-readout calorimeter + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,eta:-8,phi:9,x:32:-7,y:-7,c:1,module:2 + + + + + + + diff --git a/Detector/DRcalo/compact/DRcalo_IDEA.xml b/Detector/DRcalo/compact/DRcalo_IDEA.xml new file mode 100644 index 00000000..d80b4085 --- /dev/null +++ b/Detector/DRcalo/compact/DRcalo_IDEA.xml @@ -0,0 +1,729 @@ + + + + + + + + + + + + The compact format for the dual-readout calorimeter + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,eta:-8,phi:9,x:32:-11,y:-9,c:1,module:2 + + + + + + + diff --git a/Detector/DRcalo/compact/DRcalo_Wratten9_S13615-1025.xml b/Detector/DRcalo/compact/DRcalo_Wratten9_S13615-1025.xml new file mode 100644 index 00000000..8a145016 --- /dev/null +++ b/Detector/DRcalo/compact/DRcalo_Wratten9_S13615-1025.xml @@ -0,0 +1,740 @@ + + + + + + + + + + + + The compact format for the dual-readout calorimeter + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,eta:-8,phi:9,x:32:-7,y:-7,c:1,module:2 + + + + + + + diff --git a/Detector/DRcalo/include/DRconstructor.h b/Detector/DRcalo/include/DRconstructor.h new file mode 100644 index 00000000..dc55b446 --- /dev/null +++ b/Detector/DRcalo/include/DRconstructor.h @@ -0,0 +1,77 @@ +#ifndef DRconstructor_h +#define DRconstructor_h 1 + +#include "DRparamBarrel.h" +#include "GridDRcaloHandle.h" + +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/OpticalSurfaces.h" +#include "DD4hep/Printout.h" +#include "DD4hep/Detector.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* paramBarrel) { fParamBarrel = paramBarrel; } + void setDRparamEndcap(dd4hep::DDSegmentation::DRparamEndcap* 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 setMirrorSurf(dd4hep::OpticalSurface* mirrorSurf) { fMirrorSurf = mirrorSurf; } + void setSensDet(dd4hep::SensitiveDetector* sensDet) { + fSensDet = sensDet; + fSegmentation = dynamic_cast( sensDet->readout().segmentation().segmentation() ); + } + + void construct(); + + private: + void implementTowers(xml_comp_t& x_theta, dd4hep::DDSegmentation::DRparamBase* param); + void placeAssembly(xml_comp_t& x_theta, xml_comp_t& x_wafer, dd4hep::DDSegmentation::DRparamBase* param, + dd4hep::Trap& assemblyEnvelop, dd4hep::Volume& towerVol, dd4hep::Volume& sipmLayerVol, 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* param, int towerNo); + void implementFiber(dd4hep::Volume& towerVol, dd4hep::Trap& trap, dd4hep::Position pos, int col, int row, + dd4hep::Tube& fiberEnv, dd4hep::Tube& fiber, dd4hep::Tube& fiberC, dd4hep::Tube& fiberS, + dd4hep::Volume& capC, dd4hep::Volume& capS); + void implementSipms(dd4hep::Volume& sipmLayerVol, dd4hep::Trap& trap); + 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_hole; + xml_comp_t fX_dark; + xml_comp_t fX_mirror; + dd4hep::Assembly* fExperimentalHall; + dd4hep::Detector* fDescription; + dd4hep::DDSegmentation::DRparamBarrel* fParamBarrel; + dd4hep::DDSegmentation::DRparamEndcap* fParamEndcap; + dd4hep::DetElement* fDetElement; + dd4hep::SensitiveDetector* fSensDet; + dd4hep::OpticalSurface* fSipmSurf; + dd4hep::OpticalSurface* fMirrorSurf; + dd4hep::DDSegmentation::GridDRcalo* fSegmentation; + + bool fVis; + int fNumx, fNumy; + std::vector< std::pair > fFiberCoords; + }; +} + +#endif diff --git a/Detector/DRcalo/src/DRconstructor.cpp b/Detector/DRcalo/src/DRconstructor.cpp new file mode 100644 index 00000000..641b606e --- /dev/null +++ b/Detector/DRcalo/src/DRconstructor.cpp @@ -0,0 +1,440 @@ +#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_hole( fX_struct.child( _Unicode(hole) ) ), + fX_dark( fX_struct.child( _Unicode(dark) ) ), + fX_mirror( fX_struct.child( _Unicode(mirror) ) ) { + fExperimentalHall = nullptr; + fParamBarrel = nullptr; + fDescription = nullptr; + fDetElement = nullptr; + fSensDet = nullptr; + fSipmSurf = nullptr; + fMirrorSurf = nullptr; + fSegmentation = nullptr; + fVis = false; + fNumx = 0; + fNumy = 0; + fFiberCoords.reserve(100000); +} + +void ddDRcalo::DRconstructor::construct() { + // set vis on/off + fVis = fDescription->visAttributes(fX_det.visStr()).showDaughters(); + + implementTowers(fX_barrel, fParamBarrel); + implementTowers(fX_endcap, fParamEndcap); +} + +void ddDRcalo::DRconstructor::implementTowers(xml_comp_t& x_theta, dd4hep::DDSegmentation::DRparamBase* param) { + 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 assemblyEnvelop( (x_theta.height()+param->GetSipmHeight())/2., 0., 0., param->GetH1(), param->GetBl1(), param->GetTl1(), 0., + param->GetH2sipm(), param->GetBl2sipm(), param->GetTl2sipm(), 0. ); + + 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, towerNo); + + xml_comp_t x_wafer ( fX_sipmDim.child( _Unicode(sipmWafer) ) ); + + // Assume the top surface is nearly rectangular shape + dd4hep::Trap sipmLayer( (param->GetSipmHeight()-x_wafer.height())/2., 0., 0., param->GetH2(), param->GetBl2(), param->GetTl2(), 0., + param->GetH2(), param->GetBl2(), param->GetTl2(), 0. ); + dd4hep::Volume sipmLayerVol( "sipmLayer", sipmLayer, fDescription->material(fX_sipmDim.materialStr()) ); + if (fVis) sipmLayerVol.setVisAttributes(*fDescription, fX_sipmDim.visStr()); + + // 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); + } + + implementSipms(sipmLayerVol, tower); + + for (int nPhi = 0; nPhi < x_theta.nphi(); nPhi++) { + placeAssembly(x_theta,x_wafer,param,assemblyEnvelop,towerVol,sipmLayerVol,sipmWaferVol,towerNo,nPhi); + + if ( fX_det.reflect() ) + placeAssembly(x_theta,x_wafer,param,assemblyEnvelop,towerVol,sipmLayerVol,sipmWaferVol,towerNo,nPhi,false); + } + } + + param->filled(); + param->SetTotTowerNum( towerNo - x_theta.start() ); +} + +void ddDRcalo::DRconstructor::placeAssembly(xml_comp_t& x_theta, xml_comp_t& x_wafer, dd4hep::DDSegmentation::DRparamBase* param, + dd4hep::Trap& assemblyEnvelop, dd4hep::Volume& towerVol, dd4hep::Volume& sipmLayerVol, 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); + + // copy number of assemblyVolume is unpredictable, use dummy volume to make use of copy number of afterwards + dd4hep::Volume assemblyEnvelopVol( std::string("assembly") + (isRHS ? "" : "_refl") , assemblyEnvelop, fDescription->material("Vacuum") ); + fExperimentalHall->placeVolume( assemblyEnvelopVol, param->GetAssembleTransform3D(nPhi) ); + + assemblyEnvelopVol.placeVolume( towerVol, towerId32, dd4hep::Position(0.,0.,-param->GetSipmHeight()/2.) ); + + assemblyEnvelopVol.placeVolume( sipmLayerVol, towerId32, dd4hep::Position(0.,0.,(x_theta.height()-x_wafer.height())/2.) ); + + dd4hep::PlacedVolume sipmWaferPhys = assemblyEnvelopVol.placeVolume( sipmWaferVol, towerId32, dd4hep::Position(0.,0.,(x_theta.height()+param->GetSipmHeight()-x_wafer.height())/2.) ); + 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* param, int towerNo) { + dd4hep::Tube fiberEnv = dd4hep::Tube(0.,fX_cladC.rmax(),x_theta.height()/2.); + dd4hep::Tube fiber = dd4hep::Tube(0.,fX_cladC.rmax(),x_theta.height()/2.-fX_mirror.height()/2.); + dd4hep::Tube fiberC = dd4hep::Tube(0.,fX_coreC.rmin(),x_theta.height()/2.-fX_mirror.height()/2.); + dd4hep::Tube fiberS = dd4hep::Tube(0.,fX_coreS.rmin(),x_theta.height()/2.-fX_mirror.height()/2.); + + dd4hep::Tube cap = dd4hep::Tube(0.,fX_coreC.rmax(),fX_mirror.height()/2.); + dd4hep::Volume capC = dd4hep::Volume("capC", cap, fDescription->material(fX_mirror.materialStr())); + dd4hep::Volume capS = dd4hep::Volume("capS", cap, fDescription->material(fX_dark.materialStr())); + dd4hep::SkinSurface(*fDescription, *fDetElement, "MirrorSurf_Tower"+std::to_string(towerNo), *fMirrorSurf, capC); + if (fVis) capC.setVisAttributes(*fDescription, fX_mirror.visStr()); + if (fVis) capS.setVisAttributes(*fDescription, fX_dark.visStr()); + + 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 )/gridSize ) ) + 1; // in phi direction + fNumy = static_cast( std::floor( ( param->GetH2()*2. - sipmSize )/gridSize ) ) + 1; // in eta direction + int numxBl2 = static_cast( std::floor( ( param->GetBl2()*2. - sipmSize )/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()); + + implementFiber(unitBoxVol, trap, dd4hep::Position(-gridSize/2.,-gridSize/2.,0.), cmin, rmin, fiberEnv, fiber, fiberC, fiberS, capC, capS); + implementFiber(unitBoxVol, trap, dd4hep::Position(gridSize/2.,-gridSize/2.,0.), cmin+1, rmin, fiberEnv, fiber, fiberC, fiberS, capC, capS); + implementFiber(unitBoxVol, trap, dd4hep::Position(-gridSize/2.,gridSize/2.,0.), cmin, rmin+1, fiberEnv, fiber, fiberC, fiberS, capC, capS); + implementFiber(unitBoxVol, trap, dd4hep::Position(gridSize/2.,gridSize/2.,0.), cmin+1, rmin+1, fiberEnv, fiber, fiberC, fiberS, capC, capS); + + 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, trap, pos, column, row, fiberEnv, fiber, fiberC, fiberS, capC, capS); + 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 + checkContained(rootTrap,pos,towerHeight/2.-fiberLen,true); + + dd4hep::Position centerPos( pos.x(),pos.y(),centerZ ); + + dd4hep::Tube shortFiberEnv = dd4hep::Tube(0.,fX_cladC.rmax(),fiberLen/2.); + dd4hep::Tube shortFiber = dd4hep::Tube(0.,fX_cladC.rmax(),fiberLen/2.-fX_mirror.height()/2.); + dd4hep::Tube shortFiberC = dd4hep::Tube(0.,fX_coreC.rmin(),fiberLen/2.-fX_mirror.height()/2.); + dd4hep::Tube shortFiberS = dd4hep::Tube(0.,fX_coreS.rmin(),fiberLen/2.-fX_mirror.height()/2.); + + implementFiber(towerVol, trap, centerPos, column, row, shortFiberEnv, shortFiber, shortFiberC, shortFiberS, capC, capS); + fFiberCoords.push_back( std::make_pair(column,row) ); + } + } + } +} + +void ddDRcalo::DRconstructor::implementFiber(dd4hep::Volume& towerVol, dd4hep::Trap& trap, dd4hep::Position pos, int col, int row, + dd4hep::Tube& fiberEnv, dd4hep::Tube& fiber, dd4hep::Tube& fiberC, dd4hep::Tube& fiberS, + dd4hep::Volume& capC, dd4hep::Volume& capS) { + // punch air hole + if ( fX_hole.gap() && pos.z() > TGeoShape::Tolerance() ) { + dd4hep::Tube airHoleTube = dd4hep::Tube(0.,fX_cladC.rmax(),pos.z()); + dd4hep::Position airPos( pos.x(), pos.y(), -fiberEnv.dZ() ); + dd4hep::IntersectionSolid airHole = dd4hep::IntersectionSolid(trap,airHoleTube,airPos); + dd4hep::Volume airHoleVol("airHole", airHole, fDescription->material(fX_hole.materialStr())); + towerVol.placeVolume(airHoleVol); + } + + dd4hep::Volume fiberEnvVol("fiberEnv", fiberEnv, fDescription->material(fX_hole.materialStr())); + towerVol.placeVolume( fiberEnvVol, pos ); + + if ( fSegmentation->IsCerenkov(col,row) ) { //c fiber + dd4hep::Volume cladVol("cladC", fiber, fDescription->material(fX_cladC.materialStr())); + fiberEnvVol.placeVolume( cladVol, dd4hep::Position(0.,0.,fX_mirror.height()/2.) ); + if (fVis) cladVol.setVisAttributes(*fDescription, fX_cladC.visStr()); // high CPU consumption! + + dd4hep::Volume coreVol("coreC", fiberC, fDescription->material(fX_coreC.materialStr())); + if (fVis) coreVol.setVisAttributes(*fDescription, fX_coreC.visStr()); + cladVol.placeVolume( coreVol ); + fiberEnvVol.placeVolume( capC, dd4hep::Position(0.,0.,fX_mirror.height()/2.-fiberEnv.dZ()) ); + + coreVol.setRegion(*fDescription, fX_det.regionStr()); + cladVol.setRegion(*fDescription, fX_det.regionStr()); + } else { // s fiber + dd4hep::Volume cladVol("cladS", fiber, fDescription->material(fX_coreC.materialStr())); + fiberEnvVol.placeVolume( cladVol, dd4hep::Position(0.,0.,fX_mirror.height()/2.) ); + if (fVis) cladVol.setVisAttributes(*fDescription, fX_coreC.visStr()); + + dd4hep::Volume coreVol("coreS", fiberS, fDescription->material(fX_coreS.materialStr())); + if (fVis) coreVol.setVisAttributes(*fDescription, fX_coreS.visStr()); + cladVol.placeVolume( coreVol ); + fiberEnvVol.placeVolume( capS, dd4hep::Position(0.,0.,fX_mirror.height()/2.-fiberEnv.dZ()) ); + + coreVol.setRegion(*fDescription, fX_det.regionStr()); + cladVol.setRegion(*fDescription, fX_det.regionStr()); + } +} + +void ddDRcalo::DRconstructor::implementSipms(dd4hep::Volume& sipmLayerVol, dd4hep::Trap& trap) { + xml_comp_t x_glass ( fX_sipmDim.child( _Unicode(sipmGlass) ) ); + xml_comp_t x_wafer ( fX_sipmDim.child( _Unicode(sipmWafer) ) ); + + float sipmSize = fX_dim.dx(); + double windowHeight = fX_sipmDim.height() - x_wafer.height(); + + int rmin = 0, rmax = 0, cmin = 0, cmax = 0; + auto rootTrap = trap.access(); + dd4hep::Box sipmFullBox = calculateFullBox(rootTrap,rmin,rmax,cmin,cmax,windowHeight/2.); + dd4hep::Volume sipmFullBoxVol("sipmFullBox",sipmFullBox,fDescription->material(fX_sipmDim.materialStr())); + + float gridSize = fX_dim.distance(); + dd4hep::Box sipmUnitBox = dd4hep::Box(gridSize,gridSize,windowHeight/2.); + dd4hep::Volume sipmUnitBoxVol("sipmUnitBox",sipmUnitBox,fDescription->material(fX_sipmDim.materialStr())); + + // Glass box + dd4hep::Box sipmEnvelop(sipmSize/2., sipmSize/2., windowHeight/2.); + dd4hep::Volume sipmEnvelopVol( "sipmEnvelop", sipmEnvelop, fDescription->material(x_glass.materialStr()) ); + + if (fVis) { + sipmFullBoxVol.setVisAttributes(*fDescription, fX_sipmDim.visStr()); + sipmUnitBoxVol.setVisAttributes(*fDescription, fX_sipmDim.visStr()); + sipmEnvelopVol.setVisAttributes(*fDescription, fX_sipmDim.visStr()); + } + + sipmUnitBoxVol.placeVolume( sipmEnvelopVol, dd4hep::Position(-gridSize/2.,-gridSize/2.,0.) ); + sipmUnitBoxVol.placeVolume( sipmEnvelopVol, dd4hep::Position(gridSize/2.,-gridSize/2.,0.) ); + sipmUnitBoxVol.placeVolume( sipmEnvelopVol, dd4hep::Position(-gridSize/2.,gridSize/2.,0.) ); + sipmUnitBoxVol.placeVolume( sipmEnvelopVol, dd4hep::Position(gridSize/2.,gridSize/2.,0.) ); + + bool isEvenRow = false, isEvenCol = false; + placeUnitBox(sipmFullBoxVol,sipmUnitBoxVol,rmin,rmax,cmin,cmax,isEvenRow,isEvenCol); + sipmLayerVol.placeVolume(sipmFullBoxVol); + + for (unsigned iFiber = 0; iFiber < fFiberCoords.size(); iFiber++) { + int column = fFiberCoords.at(iFiber).first; + int row = fFiberCoords.at(iFiber).second; + auto localPosition = fSegmentation->localPosition(fNumx,fNumy,column,row); + + if ( row >= rmin && row <= rmax && column >= cmin && column <= cmax ) { + if ( ( !isEvenRow && row==rmax ) || ( !isEvenCol && column==cmax ) ) { + dd4hep::Position pos = dd4hep::Position(localPosition.x(),localPosition.y(),0.); + sipmFullBoxVol.placeVolume( sipmEnvelopVol, pos ); + } + } else { + dd4hep::Position pos = dd4hep::Position(localPosition.x(),localPosition.y(),0.); + sipmLayerVol.placeVolume( sipmEnvelopVol, pos ); + } + } + + // clear fiber coordinate vector + fFiberCoords.clear(); +} + +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., 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; + } + } + + // assume phi symmetry + float xmax = -xmin; + cmax = fNumx-1 - cmin; + + // verify assumptions + if ( std::abs(ymax+ymin) > TGeoShape::Tolerance() ) + throw std::runtime_error("Envelop of full length fibers (fullBox) is not located at the centre of the tower!"); + + 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/DRcalo/src/DRgeo.cpp b/Detector/DRcalo/src/DRgeo.cpp new file mode 100644 index 00000000..60f68c6e --- /dev/null +++ b/Detector/DRcalo/src/DRgeo.cpp @@ -0,0 +1,74 @@ +#include "DRparamBarrel.h" +#include "DRparamEndcap.h" +#include "DRconstructor.h" +#include "GridDRcaloHandle.h" + +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/OpticalSurfaces.h" +#include "DD4hep/Printout.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"); + dd4hep::OpticalSurface mirrorSurfProp = surfMgr.opticalSurface("/world/"+name+"#MirrorSurf"); + 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.setMirrorSurf(&mirrorSurfProp); + 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(ddDRcalo, ddDRcalo::create_detector) // factory method diff --git a/Detector/DRcomponents/CMakeLists.txt b/Detector/DRcomponents/CMakeLists.txt new file mode 100644 index 00000000..501c74c6 --- /dev/null +++ b/Detector/DRcomponents/CMakeLists.txt @@ -0,0 +1,31 @@ +project(DRcomponents) + +include_directories( + ${CMAKE_CURRENT_SOURCE_DIR}/include +) + +file(GLOB sources + ${PROJECT_SOURCE_DIR}/src/*.cpp +) + +file(GLOB headers + ${PROJECT_SOURCE_DIR}/include/*.h +) + +add_library(DRcomponents SHARED ${sources} ${headers}) + +target_include_directories(DRcomponents PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${DD4hep_DIR}/include +) + +target_link_libraries( + DRcomponents + DD4hep::DDCore + DD4hep::DDG4 +) + +install(TARGETS DRcomponents + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" COMPONENT dev +) diff --git a/Detector/DRcomponents/include/GeoConstruction.h b/Detector/DRcomponents/include/GeoConstruction.h new file mode 100644 index 00000000..b83d4c2c --- /dev/null +++ b/Detector/DRcomponents/include/GeoConstruction.h @@ -0,0 +1,41 @@ +#ifndef GEOCONSTRUCTION_H +#define GEOCONSTRUCTION_H 1 + +// DD4hep +#include "DDG4/Geant4GeometryInfo.h" + +// Geant4 +#include "G4VUserDetectorConstruction.hh" + +namespace dd4hep { +class Detector; +} +/** @class GeoConstruction DetectorDescription/DetDesServices/src/GeoConstruction.h GeoConstruction.h + * + * Class to create Geant4 detector geometry from TGeo representation + * On demand (ie. when calling "Construct") the DD4hep geometry is converted + * to Geant4 with all volumes, assemblies, shapes, materials etc. + * + * @author Markus Frank + * @author Anna Zaborowska + */ + +namespace det { +class GeoConstruction : public G4VUserDetectorConstruction { +public: + /// Constructor + GeoConstruction(dd4hep::Detector& lcdd); + /// Default destructor + virtual ~GeoConstruction(); + /// Geometry construction callback: Invoke the conversion to Geant4 + /// All volumes (including world) are deleted in ~G4PhysicalVolumeStore() + virtual G4VPhysicalVolume* Construct() final; + /// Construct SD + virtual void ConstructSDandField() final; + +private: + /// Reference to geometry object + dd4hep::Detector& m_lcdd; +}; +} +#endif diff --git a/Detector/DRcomponents/include/GeoSvc.h b/Detector/DRcomponents/include/GeoSvc.h new file mode 100644 index 00000000..db05e523 --- /dev/null +++ b/Detector/DRcomponents/include/GeoSvc.h @@ -0,0 +1,44 @@ +#ifndef GEOSVC_H +#define GEOSVC_H 1 + +// DD4Hep +#include "DD4hep/Detector.h" + +// Geant4 +#include "G4RunManager.hh" +#include "G4VUserDetectorConstruction.hh" + +class GeoSvc { + +public: + /// Default constructor + GeoSvc(std::vector names); + + /// Destructor + virtual ~GeoSvc(); + /// Initialize function + void initialize(); + /// This function generates the DD4hep geometry + void buildDD4HepGeo(); + /// This function generates the Geant4 geometry + void buildGeant4Geo(); + // receive DD4hep Geometry + dd4hep::DetElement getDD4HepGeo(); + dd4hep::Detector* lcdd(); + // receive Geant4 Geometry + G4VUserDetectorConstruction* getGeant4Geo(); + + static GeoSvc* GetInstance(); + +private: + /// Pointer to the interface to the DD4hep geometry + dd4hep::Detector* m_dd4hepgeo; + /// Pointer to the detector construction of DDG4 + std::shared_ptr m_geant4geo; + /// XML-files with the detector description + std::vector m_xmlFileNames; + + static GeoSvc* fInstance; +}; + +#endif // GEOSVC_H diff --git a/Detector/DRcomponents/src/GeoConstruction.cpp b/Detector/DRcomponents/src/GeoConstruction.cpp new file mode 100644 index 00000000..4cbae0ab --- /dev/null +++ b/Detector/DRcomponents/src/GeoConstruction.cpp @@ -0,0 +1,80 @@ +#include "GeoConstruction.h" + +#include + +// DD4hep +#include "DD4hep/Detector.h" +#include "DD4hep/Plugins.h" +#include "DDG4/Geant4Converter.h" +#include "TGeoManager.h" + +// Geant4 +#include "G4PVPlacement.hh" +#include "G4SDManager.hh" +#include "G4VSensitiveDetector.hh" + +namespace det { + +GeoConstruction::GeoConstruction(dd4hep::Detector& lcdd) : m_lcdd(lcdd) {} + +GeoConstruction::~GeoConstruction() {} + +// method borrowed from dd4hep::sim::Geant4DetectorSensitivesConstruction +// ::constructSensitives(Geant4DetectorConstructionContext* ctxt) +void GeoConstruction::ConstructSDandField() { + typedef std::set VolSet; + typedef std::map _SV; + dd4hep::sim::Geant4GeometryInfo* p = dd4hep::sim::Geant4Mapping::instance().ptr(); + _SV& vols = p->sensitives; + + for (_SV::const_iterator iv = vols.begin(); iv != vols.end(); ++iv) { + dd4hep::SensitiveDetector sd = (*iv).first; + std::string typ = sd.type(), nam = sd.name(); + // Sensitive detectors are deleted in ~G4SDManager + G4VSensitiveDetector* g4sd = dd4hep::PluginService::Create(typ, nam, &m_lcdd); + if (g4sd == nullptr) { + std::string tmp = typ; + tmp[0] = ::toupper(tmp[0]); + typ = "Geant4" + tmp; + g4sd = dd4hep::PluginService::Create(typ, nam, &m_lcdd); + if (g4sd == nullptr) { + dd4hep::PluginDebug dbg; + g4sd = dd4hep::PluginService::Create(typ, nam, &m_lcdd); + if (g4sd == nullptr) { + throw std::runtime_error("ConstructSDandField: FATAL Failed to " + "create Geant4 sensitive detector " + + nam + " of type " + typ + "."); + } + } + } + g4sd->Activate(true); + G4SDManager::GetSDMpointer()->AddNewDetector(g4sd); + const VolSet& sens_vols = (*iv).second; + for (VolSet::const_iterator i = sens_vols.begin(); i != sens_vols.end(); ++i) { + const TGeoVolume* vol = *i; + G4LogicalVolume* g4v = p->g4Volumes[vol]; + if (g4v == nullptr) { + throw std::runtime_error("ConstructSDandField: Failed to access G4LogicalVolume for SD " + nam + " of type " + + typ + "."); + } + G4SDManager::GetSDMpointer()->AddNewDetector(g4sd); + g4v->SetSensitiveDetector(g4sd); + } + } +} + +// method borrowed from dd4hep::sim::Geant4DetectorConstruction::Construct() +G4VPhysicalVolume* GeoConstruction::Construct() { + dd4hep::sim::Geant4Mapping& g4map = dd4hep::sim::Geant4Mapping::instance(); + dd4hep::DetElement world = m_lcdd.world(); + dd4hep::sim::Geant4Converter conv(m_lcdd, dd4hep::DEBUG); + dd4hep::sim::Geant4GeometryInfo* geo_info = conv.create(world).detach(); + g4map.attach(geo_info); + // All volumes are deleted in ~G4PhysicalVolumeStore() + G4VPhysicalVolume* m_world = geo_info->world(); + m_lcdd.apply("DD4hepVolumeManager", 0, 0); + // Create Geant4 volume manager + g4map.volumeManager(); + return m_world; +} +} diff --git a/Detector/DRcomponents/src/GeoSvc.cpp b/Detector/DRcomponents/src/GeoSvc.cpp new file mode 100644 index 00000000..ffb75bf7 --- /dev/null +++ b/Detector/DRcomponents/src/GeoSvc.cpp @@ -0,0 +1,69 @@ +#include "GeoSvc.h" +#include "GeoConstruction.h" +#include "TGeoManager.h" + +#include + +#include "DD4hep/Printout.h" + +GeoSvc* GeoSvc::fInstance = 0; + +GeoSvc::GeoSvc(std::vector names) +: m_dd4hepgeo(0), m_geant4geo(0), m_xmlFileNames(names) { + initialize(); + + if (fInstance==0) { + fInstance = this; + } +} + +GeoSvc::~GeoSvc() { + if (m_dd4hepgeo){ + try { + m_dd4hepgeo->destroyInstance(); + m_dd4hepgeo = 0; + } catch(...) {} + } +} + +void GeoSvc::initialize() { + buildDD4HepGeo(); + buildGeant4Geo(); + + return; +} + +GeoSvc* GeoSvc::GetInstance() { + return fInstance; +} + +void GeoSvc::buildDD4HepGeo() { + // we retrieve the the static instance of the DD4HEP::Geometry + m_dd4hepgeo = &(dd4hep::Detector::getInstance()); + + // load geometry + if (m_xmlFileNames.size()==0) throw std::runtime_error("List of xml file names is empty!"); + + for (auto& filename : m_xmlFileNames) { + std::cout << "loading geometry from file: '" << filename << "'" << std::endl; + m_dd4hepgeo->fromCompact(filename); + } + m_dd4hepgeo->volumeManager(); + m_dd4hepgeo->apply("DD4hepVolumeManager", 0, 0); + + return; +} + +dd4hep::Detector* GeoSvc::lcdd() { return (m_dd4hepgeo); } + +dd4hep::DetElement GeoSvc::getDD4HepGeo() { return (lcdd()->world()); } + +void GeoSvc::buildGeant4Geo() { + std::shared_ptr detector(new det::GeoConstruction(*lcdd())); + m_geant4geo = detector; + + if (m_geant4geo) return; + else throw std::runtime_error("Failed to build Geant4Geo"); +} + +G4VUserDetectorConstruction* GeoSvc::getGeant4Geo() { return (m_geant4geo.get()); } diff --git a/Detector/DRsegmentation/CMakeLists.txt b/Detector/DRsegmentation/CMakeLists.txt new file mode 100644 index 00000000..2384886e --- /dev/null +++ b/Detector/DRsegmentation/CMakeLists.txt @@ -0,0 +1,30 @@ +project(DRsegmentation) + +file(GLOB sources + ${PROJECT_SOURCE_DIR}/src/*.cpp +) + +file(GLOB headers + ${PROJECT_SOURCE_DIR}/include/*.h +) + +add_library(DRsegmentation SHARED ${sources} ${headers}) + +target_include_directories(DRsegmentation PUBLIC + $ + $ +) + +set_target_properties(DRsegmentation PROPERTIES PUBLIC_HEADER "${headers}") + +target_link_libraries( + DRsegmentation + DD4hep::DDCore +) + +dd4hep_generate_rootmap(DRsegmentation) + +install(TARGETS DRsegmentation EXPORT DetectorTargets + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" COMPONENT dev +) diff --git a/Detector/DRsegmentation/include/DRparamBarrel.h b/Detector/DRsegmentation/include/DRparamBarrel.h new file mode 100644 index 00000000..a3cdccfe --- /dev/null +++ b/Detector/DRsegmentation/include/DRparamBarrel.h @@ -0,0 +1,27 @@ +#ifndef DRparamBarrel_h +#define DRparamBarrel_h 1 + +#include "DRparamBase.h" + +#include "TVector3.h" +#include "DD4hep/DetFactoryHelper.h" + +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + class DRparamBarrel : public DRparamBase { + public: + DRparamBarrel(); + virtual ~DRparamBarrel(); + + virtual void SetDeltaThetaByTowerNo(int signedTowerNo, int) override; + virtual void SetThetaOfCenterByTowerNo(int signedTowerNo, int) override; + + virtual void init() override; + }; +} +} + +#endif diff --git a/Detector/DRsegmentation/include/DRparamBase.h b/Detector/DRsegmentation/include/DRparamBase.h new file mode 100644 index 00000000..f90ad513 --- /dev/null +++ b/Detector/DRsegmentation/include/DRparamBase.h @@ -0,0 +1,99 @@ +#ifndef DRparamBase_h +#define DRparamBase_h 1 + +#include "TVector3.h" +#include "DD4hep/DetFactoryHelper.h" + +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + class DRparamBase { + public: + DRparamBase(); + virtual ~DRparamBase(); + + 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/Detector/DRsegmentation/include/DRparamEndcap.h b/Detector/DRsegmentation/include/DRparamEndcap.h new file mode 100644 index 00000000..bcb8fd39 --- /dev/null +++ b/Detector/DRsegmentation/include/DRparamEndcap.h @@ -0,0 +1,27 @@ +#ifndef DRparamEndcap_h +#define DRparamEndcap_h 1 + +#include "DRparamBase.h" + +#include "TVector3.h" +#include "DD4hep/DetFactoryHelper.h" + +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + class DRparamEndcap : public DRparamBase { + public: + DRparamEndcap(); + virtual ~DRparamEndcap(); + + virtual void SetDeltaThetaByTowerNo(int signedTowerNo, int BEtrans) override; + virtual void SetThetaOfCenterByTowerNo(int signedTowerNo, int BEtrans) override; + + virtual void init() override; + }; +} +} + +#endif diff --git a/Detector/DRsegmentation/include/GridDRcalo.h b/Detector/DRsegmentation/include/GridDRcalo.h new file mode 100644 index 00000000..60dd787d --- /dev/null +++ b/Detector/DRsegmentation/include/GridDRcalo.h @@ -0,0 +1,109 @@ +#ifndef GridDRcalo_h +#define GridDRcalo_h 1 + +#include "DRparamBarrel.h" +#include "DRparamEndcap.h" + +#include "DDSegmentation/Segmentation.h" + +namespace dd4hep { +namespace DDSegmentation { +class GridDRcalo : public Segmentation { +public: + /// default constructor using an arbitrary type + GridDRcalo(const std::string& aCellEncoding); + /// Default constructor used by derived classes passing an existing decoder + GridDRcalo(const BitFieldCoder* decoder); + /// destructor + virtual ~GridDRcalo() override; + + // Determine the global(local) position based on the cell ID. + virtual Vector3D position(const CellID& aCellID) const; + 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; + + 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& 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 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* paramBarrel() { return fParamBarrel; } + DRparamEndcap* paramEndcap() { return fParamEndcap; } + + DRparamBase* setParamBase(int noEta) const; + +protected: + std::string fNumEtaId; + std::string fNumPhiId; + std::string fXId; + std::string fYId; + std::string fIsCerenkovId; + std::string fModule; + + double fGridSize; + double fSipmSize; + +private: + DRparamBarrel* fParamBarrel; + DRparamEndcap* fParamEndcap; +}; +} +} + +#endif diff --git a/Detector/DRsegmentation/include/GridDRcaloHandle.h b/Detector/DRsegmentation/include/GridDRcaloHandle.h new file mode 100644 index 00000000..0fb81806 --- /dev/null +++ b/Detector/DRsegmentation/include/GridDRcaloHandle.h @@ -0,0 +1,90 @@ +#ifndef GridDRcaloHandle_h +#define GridDRcaloHandle_h 1 + +#include "GridDRcalo.h" + +#include "DD4hep/Segmentations.h" +#include "DD4hep/detail/SegmentationsInterna.h" + +namespace dd4hep { + class Segmentation; + template class SegmentationWrapper; + + typedef Handle> GridDRcaloHandle; + + class GridDRcalo : public GridDRcaloHandle { + typedef GridDRcaloHandle::Object Object; + + public: + /// Default constructor + GridDRcalo() = default; + /// Copy constructor + GridDRcalo(const GridDRcalo& e) = default; + /// Copy Constructor from segmentation base object + GridDRcalo(const Segmentation& e) : Handle(e) {} + /// Copy constructor from handle + GridDRcalo(const Handle& e) : Handle(e) {} + /// Copy constructor from other polymorph/equivalent handle + template + GridDRcalo(const Handle& e) : Handle(e) {} + /// Assignment operator + GridDRcalo& operator=(const GridDRcalo& seg) = default; + /// Equality operator + bool operator==(const GridDRcalo& 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/Detector/DRsegmentation/src/DRparamBarrel.cpp b/Detector/DRsegmentation/src/DRparamBarrel.cpp new file mode 100644 index 00000000..8bce2140 --- /dev/null +++ b/Detector/DRsegmentation/src/DRparamBarrel.cpp @@ -0,0 +1,98 @@ +#include "DRparamBarrel.h" + +#include "Math/GenVector/RotationZYX.h" + +#include + +namespace dd4hep { +namespace DDSegmentation { + +DRparamBarrel::DRparamBarrel() { + 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::~DRparamBarrel() {} + +void DRparamBarrel::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::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::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/Detector/DRsegmentation/src/DRparamBase.cpp b/Detector/DRsegmentation/src/DRparamBase.cpp new file mode 100644 index 00000000..5630449e --- /dev/null +++ b/Detector/DRsegmentation/src/DRparamBase.cpp @@ -0,0 +1,100 @@ +#include "DRparamBase.h" + +#include "Math/GenVector/RotationZYX.h" + +#include + +namespace dd4hep { +namespace DDSegmentation { + +DRparamBase::DRparamBase() { + 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::~DRparamBase() {} + +dd4hep::RotationZYX DRparamBase::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::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::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::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::GetTransform3D(int numPhi) { + auto rot = GetRotationZYX(numPhi); + auto pos = GetTowerPos(numPhi); + + return dd4hep::Transform3D(rot,pos); +} + +dd4hep::Transform3D DRparamBase::GetAssembleTransform3D(int numPhi) { + auto rot = GetRotationZYX(numPhi); + auto pos = GetAssemblePos(numPhi); + + return dd4hep::Transform3D(rot,pos); +} + +dd4hep::Transform3D DRparamBase::GetSipmTransform3D(int numPhi) { + auto rot = GetRotationZYX(numPhi); + auto pos = GetSipmLayerPos(numPhi); + + return dd4hep::Transform3D(rot,pos); +} + +} +} diff --git a/Detector/DRsegmentation/src/DRparamEndcap.cpp b/Detector/DRsegmentation/src/DRparamEndcap.cpp new file mode 100644 index 00000000..9e554f7f --- /dev/null +++ b/Detector/DRsegmentation/src/DRparamEndcap.cpp @@ -0,0 +1,98 @@ +#include "DRparamEndcap.h" + +#include "Math/GenVector/RotationZYX.h" + +#include + +namespace dd4hep { +namespace DDSegmentation { + +DRparamEndcap::DRparamEndcap() { + 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::~DRparamEndcap() {} + +void DRparamEndcap::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::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::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/Detector/DRsegmentation/src/GridDRcalo.cpp b/Detector/DRsegmentation/src/GridDRcalo.cpp new file mode 100644 index 00000000..a556c5d6 --- /dev/null +++ b/Detector/DRsegmentation/src/GridDRcalo.cpp @@ -0,0 +1,232 @@ +#include "GridDRcalo.h" + +#include +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + +/// default constructor using an encoding string +GridDRcalo::GridDRcalo(const std::string& cellEncoding) : Segmentation(cellEncoding) { + // define type and description + _type = "GridDRcalo"; + _description = "DRcalo segmentation based on the tower / (Cerenkov or Scintillation) fiber / SiPM hierarchy"; + + // register all necessary parameters + 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(); + fParamEndcap = new DRparamEndcap(); +} + +GridDRcalo::GridDRcalo(const BitFieldCoder* decoder) : Segmentation(decoder) { + // define type and description + _type = "GridDRcalo"; + _description = "DRcalo segmentation based on the tower / (Cerenkov or Scintillation) fiber / SiPM hierarchy"; + + // register all necessary parameters + 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(); + fParamEndcap = new DRparamEndcap(); +} + +GridDRcalo::~GridDRcalo() { + if (fParamBarrel) delete fParamBarrel; + if (fParamEndcap) delete fParamEndcap; +} + +Vector3D GridDRcalo::position(const CellID& cID) const { + int noEta = numEta(cID); + int noPhi = numPhi(cID); + + DRparamBase* 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::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::localPosition(int numx, int numy, int x_, int y_) const { + 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::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; + + return setCellID( numEta(vID), numPhi(vID), x, y ); +} + +VolumeID GridDRcalo::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::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); + + return vID; +} + +// Get the identifier number of a mother tower in eta or phi direction +int GridDRcalo::numEta(const CellID& aCellID) const { + VolumeID numEta = static_cast(_decoder->get(aCellID, fNumEtaId)); + return static_cast(numEta); +} + +int GridDRcalo::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::numX(const CellID& aCellID) const { + int noEta = numEta(aCellID); + + DRparamBase* paramBase = setParamBase(noEta); + + int noX = static_cast( std::floor( ( paramBase->GetTl2()*2. - fSipmSize )/fGridSize ) ) + 1; // in phi direction + + return noX; +} + +int GridDRcalo::numY(const CellID& aCellID) const { + int noEta = numEta(aCellID); + + DRparamBase* paramBase = setParamBase(noEta); + + int noY = static_cast( std::floor( ( paramBase->GetH2()*2. - fSipmSize )/fGridSize ) ) + 1; // in eta direction + + return noY; +} + +// Get the identifier number of a SiPM in x or y direction (local coordinate) +int GridDRcalo::x(const CellID& aCellID) const { // approx eta direction + VolumeID x = static_cast(_decoder->get(aCellID, fXId)); + return static_cast(x); +} +int GridDRcalo::y(const CellID& aCellID) const { // approx phi direction + VolumeID y = static_cast(_decoder->get(aCellID, fYId)); + return static_cast(y); +} + +bool GridDRcalo::IsCerenkov(const CellID& aCellID) const { + VolumeID isCeren = static_cast(_decoder->get(aCellID, fIsCerenkovId)); + return static_cast(isCeren); +} + +bool GridDRcalo::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::IsTower(const CellID& aCellID) const { + VolumeID module = static_cast(_decoder->get(aCellID, fModule)); + return module==0; +} + +bool GridDRcalo::IsSiPM(const CellID& aCellID) const { + VolumeID module = static_cast(_decoder->get(aCellID, fModule)); + return module==1; +} + +int GridDRcalo::getLast32bits(const CellID& aCellID) const { + CellID aId64 = aCellID >> sizeof(int)*CHAR_BIT; + int aId32 = (int)aId64; + + return aId32; +} + +CellID GridDRcalo::convertLast32to64(const int aId32) const { + CellID aId64 = (CellID)aId32; + aId64 <<= sizeof(int)*CHAR_BIT; + + return aId64; +} + +DRparamBase* GridDRcalo::setParamBase(int noEta) const { + DRparamBase* 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::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/Detector/DRsegmentation/src/GridDRcaloHandle.cpp b/Detector/DRsegmentation/src/GridDRcaloHandle.cpp new file mode 100644 index 00000000..6450a673 --- /dev/null +++ b/Detector/DRsegmentation/src/GridDRcaloHandle.cpp @@ -0,0 +1,4 @@ +#include "GridDRcaloHandle.h" +#include "DD4hep/detail/Handle.inl" + +DD4HEP_INSTANTIATE_HANDLE_UNNAMED(dd4hep::DDSegmentation::GridDRcalo); diff --git a/Detector/DRsegmentation/src/SegmentationFactories.cpp b/Detector/DRsegmentation/src/SegmentationFactories.cpp new file mode 100644 index 00000000..71f4177c --- /dev/null +++ b/Detector/DRsegmentation/src/SegmentationFactories.cpp @@ -0,0 +1,12 @@ +#include "DD4hep/Factories.h" +#include "DD4hep/detail/SegmentationsInterna.h" + +namespace { +template +dd4hep::SegmentationObject* create_segmentation(const dd4hep::BitFieldCoder* decoder) { + return new dd4hep::SegmentationWrapper(decoder); +} +} + +#include "GridDRcalo.h" +DECLARE_SEGMENTATION(GridDRcalo, create_segmentation) diff --git a/Detector/DRsensitive/CMakeLists.txt b/Detector/DRsensitive/CMakeLists.txt new file mode 100644 index 00000000..4a830efe --- /dev/null +++ b/Detector/DRsensitive/CMakeLists.txt @@ -0,0 +1,32 @@ +project(DRsensitive) + +file(GLOB sources + ${PROJECT_SOURCE_DIR}/src/*.cpp +) + +file(GLOB headers + ${PROJECT_SOURCE_DIR}/include/*.h +) + +add_library(DRsensitive SHARED ${sources} ${headers}) + +target_include_directories(DRsensitive PUBLIC + $ + $ +) + +set_target_properties(DRsensitive PROPERTIES PUBLIC_HEADER "${headers}") + +target_link_libraries( + DRsensitive + DRsegmentation + DD4hep::DDCore + DD4hep::DDG4 +) + +dd4hep_generate_rootmap(DRsensitive) + +install(TARGETS DRsensitive EXPORT DetectorTargets + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" COMPONENT dev +) diff --git a/Detector/DRsensitive/include/DRcaloSiPMHit.h b/Detector/DRsensitive/include/DRcaloSiPMHit.h new file mode 100644 index 00000000..bb764d15 --- /dev/null +++ b/Detector/DRsensitive/include/DRcaloSiPMHit.h @@ -0,0 +1,68 @@ +#ifndef DRcaloSiPMHit_h +#define DRcaloSiPMHit_h 1 + +#include "G4VHit.hh" +#include "G4THitsCollection.hh" +#include "G4Allocator.hh" +#include "G4ThreeVector.hh" + +#include "DD4hep/Objects.h" +#include "DD4hep/Segmentations.h" + +namespace drc { + class DRcaloSiPMHit : public G4VHit { + public: + typedef std::map DRsimTimeStruct; + typedef std::map DRsimWavlenSpectrum; + + DRcaloSiPMHit(float wavSampling, float timeSampling); + DRcaloSiPMHit(const DRcaloSiPMHit &right); + virtual ~DRcaloSiPMHit(); + + const DRcaloSiPMHit& operator=(const DRcaloSiPMHit &right); + G4bool operator==(const DRcaloSiPMHit &right) const; + + inline void *operator new(size_t); + inline void operator delete(void* aHit); + + virtual void Draw() {}; + virtual void Print() {}; + + 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(float center); + const DRsimWavlenSpectrum& GetWavlenSpectrum() const { return fWavlenSpectrum; } + + void CountTimeStruct(float center); + const DRsimTimeStruct& GetTimeStruct() const { return fTimeStruct; } + + float GetSamplingTime() { return mTimeSampling; } + float GetSamplingWavlen() { return mWavSampling; } + + private: + dd4hep::DDSegmentation::CellID fSiPMnum; + unsigned long fPhotons; + DRsimWavlenSpectrum fWavlenSpectrum; + DRsimTimeStruct fTimeStruct; + float mWavSampling; + float mTimeSampling; + }; + + typedef G4THitsCollection DRcaloSiPMHitsCollection; + extern G4ThreadLocal G4Allocator* DRcaloSiPMHitAllocator; + + inline void* DRcaloSiPMHit::operator new(size_t) { + if (!DRcaloSiPMHitAllocator) DRcaloSiPMHitAllocator = new G4Allocator; + return (void*)DRcaloSiPMHitAllocator->MallocSingle(); + } + + inline void DRcaloSiPMHit::operator delete(void*aHit) { + DRcaloSiPMHitAllocator->FreeSingle((DRcaloSiPMHit*) aHit); + } +} + +#endif diff --git a/Detector/DRsensitive/include/DRcaloSiPMSD.h b/Detector/DRsensitive/include/DRcaloSiPMSD.h new file mode 100644 index 00000000..3e5939c4 --- /dev/null +++ b/Detector/DRsensitive/include/DRcaloSiPMSD.h @@ -0,0 +1,44 @@ +#ifndef DRcaloSiPMSD_h +#define DRcaloSiPMSD_h 1 + +#include "DRcaloSiPMHit.h" +#include "GridDRcalo.h" +#include "DD4hep/Segmentations.h" + +#include "G4VSensitiveDetector.hh" +#include "G4SystemOfUnits.hh" +#include "G4PhysicalConstants.hh" +#include "G4Step.hh" +#include "G4TouchableHistory.hh" + +namespace drc { + class DRcaloSiPMSD : public G4VSensitiveDetector { + public: + DRcaloSiPMSD(const std::string aName, const std::string aReadoutName, const dd4hep::Segmentation& aSeg); + ~DRcaloSiPMSD(); + + virtual void Initialize(G4HCofThisEvent* HCE) final; + virtual bool ProcessHits(G4Step* aStep, G4TouchableHistory*) final; + + private: + DRcaloSiPMHitsCollection* fHitCollection; + dd4hep::DDSegmentation::GridDRcalo* fSeg; + G4int fHCID; + + G4int fWavBin; + G4int fTimeBin; + G4float fWavlenStart; + G4float fWavlenEnd; + G4float fTimeStart; + G4float fTimeEnd; + G4float fWavlenStep; + G4float fTimeStep; + + G4double wavToE(G4double wav) { return h_Planck*c_light/wav; } + + float findWavCenter(G4double en); + float findTimeCenter(G4double stepTime); + }; +} + +#endif diff --git a/Detector/DRsensitive/src/DRcaloSiPMHit.cpp b/Detector/DRsensitive/src/DRcaloSiPMHit.cpp new file mode 100644 index 00000000..9ea7d76d --- /dev/null +++ b/Detector/DRsensitive/src/DRcaloSiPMHit.cpp @@ -0,0 +1,49 @@ +#include "DRcaloSiPMHit.h" + +G4ThreadLocal G4Allocator* drc::DRcaloSiPMHitAllocator = 0; + +drc::DRcaloSiPMHit::DRcaloSiPMHit(float wavSampling, float timeSampling) +: G4VHit(), + fSiPMnum(0), + fPhotons(0), + mWavSampling(wavSampling), + mTimeSampling(timeSampling) +{} + +drc::DRcaloSiPMHit::~DRcaloSiPMHit() {} + +drc::DRcaloSiPMHit::DRcaloSiPMHit(const drc::DRcaloSiPMHit &right) +: G4VHit() { + fSiPMnum = right.fSiPMnum; + fPhotons = right.fPhotons; + fWavlenSpectrum = right.fWavlenSpectrum; + fTimeStruct = right.fTimeStruct; + mWavSampling = right.mWavSampling; + mTimeSampling = right.mTimeSampling; +} + +const drc::DRcaloSiPMHit& drc::DRcaloSiPMHit::operator=(const drc::DRcaloSiPMHit &right) { + fSiPMnum = right.fSiPMnum; + fPhotons = right.fPhotons; + fWavlenSpectrum = right.fWavlenSpectrum; + fTimeStruct = right.fTimeStruct; + mWavSampling = right.mWavSampling; + mTimeSampling = right.mTimeSampling; + return *this; +} + +G4bool drc::DRcaloSiPMHit::operator==(const drc::DRcaloSiPMHit &right) const { + return (fSiPMnum==right.fSiPMnum); +} + +void drc::DRcaloSiPMHit::CountWavlenSpectrum(float center) { + auto it = fWavlenSpectrum.find(center); + if (it==fWavlenSpectrum.end()) fWavlenSpectrum.insert(std::make_pair(center,1)); + else it->second++; +} + +void drc::DRcaloSiPMHit::CountTimeStruct(float center) { + auto it = fTimeStruct.find(center); + if (it==fTimeStruct.end()) fTimeStruct.insert(std::make_pair(center,1)); + else it->second++; +} diff --git a/Detector/DRsensitive/src/DRcaloSiPMSD.cpp b/Detector/DRsensitive/src/DRcaloSiPMSD.cpp new file mode 100644 index 00000000..65252b74 --- /dev/null +++ b/Detector/DRsensitive/src/DRcaloSiPMSD.cpp @@ -0,0 +1,101 @@ +#include "DRcaloSiPMSD.h" +#include "DRcaloSiPMHit.h" + +#include "DDG4/Geant4Mapping.h" +#include "DDG4/Geant4VolumeManager.h" + +#include "G4HCofThisEvent.hh" +#include "G4SDManager.hh" +#include "G4ParticleDefinition.hh" +#include "G4ParticleTypes.hh" + +#include "G4SystemOfUnits.hh" +#include "DD4hep/DD4hepUnits.h" + +drc::DRcaloSiPMSD::DRcaloSiPMSD(const std::string aName, const std::string aReadoutName, const dd4hep::Segmentation& aSeg) +: G4VSensitiveDetector(aName), fHitCollection(0), fHCID(-1), +fWavBin(120), fTimeBin(600), fWavlenStart(900.), fWavlenEnd(300.), fTimeStart(10.), fTimeEnd(70.) +{ + collectionName.insert(aReadoutName); + fSeg = dynamic_cast( aSeg.segmentation() ); + fWavlenStep = (fWavlenStart-fWavlenEnd)/(float)fWavBin; + fTimeStep = (fTimeEnd-fTimeStart)/(float)fTimeBin; +} + +drc::DRcaloSiPMSD::~DRcaloSiPMSD() {} + +void drc::DRcaloSiPMSD::Initialize(G4HCofThisEvent* hce) { + fHitCollection = new drc::DRcaloSiPMHitsCollection(SensitiveDetectorName,collectionName[0]); + if (fHCID<0) { fHCID = GetCollectionID(0); } + hce->AddHitsCollection(fHCID,fHitCollection); +} + +G4bool drc::DRcaloSiPMSD::ProcessHits(G4Step* step, G4TouchableHistory*) { + if(step->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return false; + + 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 ); + 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 = fSeg->cellID(loc, glob, volID); + + G4int nofHits = fHitCollection->entries(); + G4double hitTime = step->GetPostStepPoint()->GetGlobalTime(); + G4double energy = step->GetTrack()->GetTotalEnergy(); + + drc::DRcaloSiPMHit* hit = NULL; + + for (G4int i = 0; i < nofHits; i++) { + if ( (*fHitCollection)[i]->GetSiPMnum()==cID ) { + hit = (*fHitCollection)[i]; + break; + } + } + + if (hit==NULL) { + hit = new DRcaloSiPMHit(fWavlenStep,fTimeStep); + hit->SetSiPMnum(cID); + + fHitCollection->insert(hit); + } + + hit->photonCount(); + + float wavCenter = findWavCenter(energy); + hit->CountWavlenSpectrum(wavCenter); + + float timeCenter = findTimeCenter(hitTime); + hit->CountTimeStruct(timeCenter); + + return true; +} + +float drc::DRcaloSiPMSD::findWavCenter(G4double en) { + int i = 0; + for ( ; i < fWavBin+1; i++) { + if ( en < wavToE( (fWavlenStart - static_cast(i)*fWavlenStep)*nm ) ) break; + } + + if (i==0) return (fWavlenStart + 0.5*fWavlenStep); + else if (i==fWavBin+1) return (fWavlenEnd - 0.5*fWavlenStep); + + return ( fWavlenStart-static_cast(i)*fWavlenStep + fWavlenStart-static_cast(i-1)*fWavlenStep )/2.; +} + +float drc::DRcaloSiPMSD::findTimeCenter(G4double stepTime) { + int i = 0; + for ( ; i < fTimeBin+1; i++) { + if ( stepTime < ( (fTimeStart + static_cast(i)*fTimeStep)*CLHEP::ns ) ) break; + } + + if (i==0) return (fTimeStart - 0.5*fTimeStep); + else if (i==fTimeBin+1) return (fTimeEnd + 0.5*fTimeStep); + + return ( fTimeStart+static_cast(i-1)*fTimeStep + fTimeStart+static_cast(i)*fTimeStep )/2.; +} diff --git a/Detector/DRsensitive/src/SDWrapper.cpp b/Detector/DRsensitive/src/SDWrapper.cpp new file mode 100644 index 00000000..fdc020ff --- /dev/null +++ b/Detector/DRsensitive/src/SDWrapper.cpp @@ -0,0 +1,14 @@ +#include "DRcaloSiPMSD.h" + +#include "DD4hep/Detector.h" +#include "DDG4/Factories.h" + +namespace dd4hep { +namespace sim { + static G4VSensitiveDetector* create_DRcaloSiPM_sd(const std::string& aDetectorName, dd4hep::Detector& aLcdd) { + std::string readoutName = aLcdd.sensitiveDetector(aDetectorName).readout().name(); + return new drc::DRcaloSiPMSD(aDetectorName,readoutName,aLcdd.sensitiveDetector(aDetectorName).readout().segmentation()); + } +} +} +DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(DRcaloSiPMSD,dd4hep::sim::create_DRcaloSiPM_sd)