diff --git a/CMakeLists.txt b/CMakeLists.txt
index a5fc00744..c02ec4a9f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -68,6 +68,7 @@ endif()
file(GLOB sources
./detector/tracker/*.cpp
./detector/calorimeter/*.cpp
+ ./detector/calorimeter/dual-readout/src/*.cpp
./detector/fcal/*.cpp
./detector/other/*.cpp
./detector/CaloTB/*.cpp
@@ -106,9 +107,15 @@ if(NOT DCH_INFO_H_EXIST)
message(WARNING "Subdetector ${FILES_DEPENDINGON_DCH_INFO_H} will not be built because header file ${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)
@@ -124,14 +131,18 @@ 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)
target_link_libraries(${PackageName}G4 LCIO::lcio)
endif()
+
#Create this_package.sh file, and install
dd4hep_instantiate_package(${PackageName})
diff --git a/FCCee/IDEA/compact/IDEA_o1_v02/DectDimensions_IDEA_o1_v01.xml b/FCCee/IDEA/compact/IDEA_o1_v02/DectDimensions_IDEA_o1_v01.xml
index 90dc7d53a..04b8a020a 100644
--- a/FCCee/IDEA/compact/IDEA_o1_v02/DectDimensions_IDEA_o1_v01.xml
+++ b/FCCee/IDEA/compact/IDEA_o1_v02/DectDimensions_IDEA_o1_v01.xml
@@ -103,45 +103,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -212,6 +173,7 @@
+
@@ -228,6 +190,8 @@
+
+
diff --git a/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml b/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml
index 4619c7cfb..7716d8087 100644
--- a/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml
+++ b/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml
@@ -52,7 +52,7 @@
-
+
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 f3779cf19..090f19d87 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,6 +53,9 @@
+
+
+
@@ -105,45 +108,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -214,6 +178,15 @@
+
+
+
+
+
+
+
+
+
@@ -230,6 +203,8 @@
+
+
@@ -262,6 +237,15 @@
+
+
+
+
+
+
+
+
+
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..8fa6606ee
--- /dev/null
+++ b/FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml
@@ -0,0 +1,703 @@
+
+
+
+
+
+
+
+ The compact format for the dual-readout calorimeter (for FCCee IDEA)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ system:5,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 ff2bbf3db..7113e8182 100644
--- a/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml
+++ b/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml
@@ -42,7 +42,7 @@
-
+
@@ -53,6 +53,9 @@
+
+
+
diff --git a/detector/calorimeter/dual-readout/include/DRconstructor.h b/detector/calorimeter/dual-readout/include/DRconstructor.h
new file mode 100644
index 000000000..eddc2b098
--- /dev/null
+++ b/detector/calorimeter/dual-readout/include/DRconstructor.h
@@ -0,0 +1,76 @@
+#ifndef DRconstructor_h
+#define DRconstructor_h 1
+
+#include "detectorSegmentations/DRparamBarrel_k4geo.h"
+#include "detectorSegmentations/GridDRcaloHandle_k4geo.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_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 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_k4geo* param);
+ void placeAssembly(xml_comp_t& x_theta, xml_comp_t& x_wafer, dd4hep::DDSegmentation::DRparamBase_k4geo* param,
+ dd4hep::Trap& assemblyEnvelop, 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, 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);
+ 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_k4geo* fParamBarrel;
+ dd4hep::DDSegmentation::DRparamEndcap_k4geo* fParamEndcap;
+ dd4hep::DetElement* fDetElement;
+ dd4hep::SensitiveDetector* fSensDet;
+ dd4hep::OpticalSurface* fSipmSurf;
+ dd4hep::OpticalSurface* fMirrorSurf;
+ dd4hep::DDSegmentation::GridDRcalo_k4geo* fSegmentation;
+
+ bool fVis;
+ int fNumx, fNumy;
+ std::vector< std::pair > fFiberCoords;
+ };
+}
+
+#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..4c6b09f56
--- /dev/null
+++ b/detector/calorimeter/dual-readout/src/DRconstructor.cpp
@@ -0,0 +1,367 @@
+#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_k4geo* 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;
+
+ std::cout << "test Tower No : " << towerNo << std::endl;
+
+ // 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) ) );
+
+ // 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(x_theta,x_wafer,param,assemblyEnvelop,towerVol,sipmWaferVol,towerNo,nPhi);
+
+ if ( fX_det.reflect() )
+ placeAssembly(x_theta,x_wafer,param,assemblyEnvelop,towerVol,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_k4geo* param,
+ dd4hep::Trap& assemblyEnvelop, 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);
+
+ // 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.) );
+
+ // Remove sipmLayer
+
+ 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_k4geo* 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.);
+
+ 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());
+
+ // // Remove cap (mirror or black paint in front of the fiber)
+ implementFiber(unitBoxVol, trap, dd4hep::Position(-gridSize/2.,-gridSize/2.,0.), cmin, rmin, fiberEnv, fiber, fiberC, fiberS);
+ implementFiber(unitBoxVol, trap, dd4hep::Position(gridSize/2.,-gridSize/2.,0.), cmin+1, rmin, fiberEnv, fiber, fiberC, fiberS);
+ implementFiber(unitBoxVol, trap, dd4hep::Position(-gridSize/2.,gridSize/2.,0.), cmin, rmin+1, fiberEnv, fiber, fiberC, fiberS);
+ implementFiber(unitBoxVol, trap, dd4hep::Position(gridSize/2.,gridSize/2.,0.), cmin+1, rmin+1, fiberEnv, fiber, fiberC, fiberS);
+
+ 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);
+ 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 );
+
+ 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);
+ 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::Trap& trap, dd4hep::Position pos, int col, int row,
+ dd4hep::Tube& fiberEnv, dd4hep::Tube& fiber, dd4hep::Tube& fiberC, dd4hep::Tube& fiberS) {
+ 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 );
+
+ 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 );
+
+ 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..8e6076265
--- /dev/null
+++ b/detector/calorimeter/dual-readout/src/FiberDualReadoutCalo_o1_v01.cpp
@@ -0,0 +1,76 @@
+#include "detectorSegmentations/DRparamBarrel_k4geo.h"
+#include "detectorSegmentations/DRparamEndcap_k4geo.h"
+#include "detectorSegmentations/GridDRcaloHandle_k4geo.h"
+
+#include "DRconstructor.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(FiberDualReadoutCalo_o1_v01, ddDRcalo::create_detector) // factory method
\ No newline at end of file
diff --git a/detectorSegmentations/include/detectorSegmentations/DRparamBarrel_k4geo.h b/detectorSegmentations/include/detectorSegmentations/DRparamBarrel_k4geo.h
new file mode 100644
index 000000000..47b2bb6af
--- /dev/null
+++ b/detectorSegmentations/include/detectorSegmentations/DRparamBarrel_k4geo.h
@@ -0,0 +1,27 @@
+#ifndef DETSEGMENTATION_DRPARAMBARREL_H
+#define DETSEGMENTATION_DRPARAMBARREL_H
+
+#include "detectorSegmentations/DRparamBase_k4geo.h"
+
+#include "TVector3.h"
+#include "DD4hep/DetFactoryHelper.h"
+
+#include
+#include
+
+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..d8fb48d4d
--- /dev/null
+++ b/detectorSegmentations/include/detectorSegmentations/DRparamEndcap_k4geo.h
@@ -0,0 +1,27 @@
+#ifndef DETSEGMENTATION_DRPARAMENDCAP_H
+#define DETSEGMENTATION_DRPARAMENDCAP_H
+
+#include "detectorSegmentations/DRparamBase_k4geo.h"
+
+#include "TVector3.h"
+#include "DD4hep/DetFactoryHelper.h"
+
+#include
+#include
+
+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/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