diff --git a/Detector/DetFCCeeECalInclined/CMakeLists.txt b/Detector/DetFCCeeECalInclined/CMakeLists.txt index 30282bcd..7282cc39 100644 --- a/Detector/DetFCCeeECalInclined/CMakeLists.txt +++ b/Detector/DetFCCeeECalInclined/CMakeLists.txt @@ -3,3 +3,9 @@ ################################################################################ install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/compact DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}/Detector/DetFCCeeECalInclined) + +file(GLOB sources src/*.cpp) +add_dd4hep_plugin(DetFCCeeECalInclined SHARED ${sources}) +target_link_libraries(DetFCCeeECalInclined DD4hep::DDCore) + +set(LIBRARY_OUTPUT_PATH ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}) diff --git a/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_upstream.xml b/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_upstream.xml index 20e4a308..b3eee7d8 100644 --- a/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_upstream.xml +++ b/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_upstream.xml @@ -11,7 +11,9 @@ - + + + @@ -132,4 +134,4 @@ - \ No newline at end of file + diff --git a/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_v2.xml b/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_v2.xml new file mode 100644 index 00000000..98ec8a50 --- /dev/null +++ b/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_v2.xml @@ -0,0 +1,163 @@ + + + + + + Settings for the inclined EM calorimeter using FCCee factory (EmCaloBarrelInclinedFCCee). + The barrel is filled with liquid argon. Passive material includes lead in the middle and steal on the outside, glued together. + Passive plates are inclined by a certain angle from the radial direction. + In between of two passive plates there is a readout. + Space between the plate and readout is of trapezoidal shape and filled with liquid argon. + Definition of sizes, visualization settings, readout and longitudinal segmentation are specified. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,cryo:1,type:4,subtype:3,layer:8,module:11,eta:9 + + + + + + + system:4,cryo:1,type:4,subtype:3,layer:8,eta:9,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Detector/DetFCCeeECalInclined/src/ECalBarrelInclinedFCCee_geo.cpp b/Detector/DetFCCeeECalInclined/src/ECalBarrelInclinedFCCee_geo.cpp new file mode 100644 index 00000000..bb224f8b --- /dev/null +++ b/Detector/DetFCCeeECalInclined/src/ECalBarrelInclinedFCCee_geo.cpp @@ -0,0 +1,628 @@ +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Handle.h" +#include "XML/Utilities.h" + +#include + +// todo: remove gaudi logging and properly capture output +#define endmsg std::endl +#define lLog std::cout +namespace MSG { +const std::string ERROR = " Error: "; +const std::string DEBUG = " Debug: "; +const std::string INFO = " Info: "; +} + +namespace det { +static dd4hep::detail::Ref_t createECalBarrelInclinedFCCee(dd4hep::Detector& aLcdd, + dd4hep::xml::Handle_t aXmlElement, + dd4hep::SensitiveDetector aSensDet) { + + dd4hep::xml::DetElement xmlDetElem = aXmlElement; + std::string nameDet = xmlDetElem.nameStr(); + dd4hep::xml::Dimension dim(xmlDetElem.dimensions()); + dd4hep::DetElement caloDetElem(nameDet, xmlDetElem.id()); + + // Create air envelope for the whole barrel + dd4hep::Volume envelopeVol(nameDet + "_vol", dd4hep::Tube(dim.rmin(), dim.rmax(), dim.dz()), + aLcdd.material("Air")); + envelopeVol.setVisAttributes(aLcdd, dim.visStr()); + + // Retrieve cryostat data + dd4hep::xml::DetElement cryostat = aXmlElement.child(_Unicode(cryostat)); + + dd4hep::xml::DetElement cryoFrontWarm = cryostat.child(_Unicode(cryo_front_warm)); + dd4hep::xml::DetElement cryoFrontCold = cryostat.child(_Unicode(cryo_front_cold)); + dd4hep::xml::DetElement cryoBackCold = cryostat.child(_Unicode(cryo_back_cold)); + dd4hep::xml::DetElement cryoSolenoid = cryostat.child(_Unicode(cryo_solenoid)); + dd4hep::xml::DetElement cryoBackWarm = cryostat.child(_Unicode(cryo_back_warm)); + dd4hep::xml::DetElement cryoSideWarm = cryostat.child(_Unicode(cryo_side_warm)); + dd4hep::xml::DetElement cryoSideCold = cryostat.child(_Unicode(cryo_side_cold)); + + bool cryoFrontSensitive = cryoFrontWarm.isSensitive(); + bool cryoBackSensitive = cryoBackWarm.isSensitive(); + bool cryoSideSensitive = cryoSideWarm.isSensitive(); + + dd4hep::xml::Dimension cryoFrontWarmDim(cryoFrontWarm.dimensions()); + // if cryoThicknessFront>0 -> build the cryostat + double cryoThicknessFront = cryoFrontWarmDim.rmax() - cryoFrontWarmDim.rmin(); + dd4hep::xml::Dimension cryoFrontColdDim(cryoFrontCold.dimensions()); + dd4hep::xml::Dimension cryoBackColdDim(cryoBackCold.dimensions()); + dd4hep::xml::Dimension cryoSolenoidDim(cryoSolenoid.dimensions()); + dd4hep::xml::Dimension cryoBackWarmDim(cryoBackWarm.dimensions()); + dd4hep::xml::Dimension cryoSideColdDim(cryoSideCold.dimensions()); + dd4hep::xml::Dimension cryoSideWarmDim(cryoSideWarm.dimensions()); + + // Retrieve active and passive material data + dd4hep::xml::DetElement calo = aXmlElement.child(_Unicode(calorimeter)); + dd4hep::xml::Dimension caloDim(calo.dimensions()); + dd4hep::xml::DetElement active = calo.child(_Unicode(active)); + std::string activeMaterial = active.materialStr(); + double activeThickness = active.thickness(); + + dd4hep::xml::DetElement overlap = active.child(_Unicode(overlap)); + double activePassiveOverlap = overlap.offset(); + if (activePassiveOverlap < 0 || activePassiveOverlap > 0.5) { + // todo: ServiceHandle incidentSvc("IncidentSvc", "ECalConstruction"); + lLog << MSG::ERROR << "Overlap between active and passive cannot be more than half of passive plane!" << endmsg; + //todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure")); + } + dd4hep::xml::DetElement layers = calo.child(_Unicode(layers)); + uint numLayers = 0; + std::vector layerHeight; + double layersTotalHeight = 0; + for (dd4hep::xml::Collection_t layer_coll(layers, _Unicode(layer)); layer_coll; ++layer_coll) { + dd4hep::xml::Component layer = layer_coll; + numLayers += layer.repeat(); + for (int iLay = 0; iLay < layer.repeat(); iLay++) { + layerHeight.push_back(layer.thickness()); + } + layersTotalHeight += layer.repeat() * layer.thickness(); + } + lLog << MSG::DEBUG << "**** FCCee ECAL Inclined ****" << endmsg; + lLog << MSG::DEBUG << "Number of layers: " << numLayers << " total thickness " << layersTotalHeight << endmsg; + + dd4hep::xml::DetElement readout = calo.child(_Unicode(readout)); + std::string readoutMaterial = readout.materialStr(); + double readoutThickness = readout.thickness(); + + dd4hep::xml::DetElement passive = calo.child(_Unicode(passive)); + dd4hep::xml::DetElement passiveInner = passive.child(_Unicode(inner)); + dd4hep::xml::DetElement passiveInnerMax = passive.child(_Unicode(innerMax)); + dd4hep::xml::DetElement passiveOuter = passive.child(_Unicode(outer)); + dd4hep::xml::DetElement passiveGlue = passive.child(_Unicode(glue)); + std::string passiveInnerMaterial = passiveInner.materialStr(); + std::string passiveOuterMaterial = passiveOuter.materialStr(); + std::string passiveGlueMaterial = passiveGlue.materialStr(); + double passiveInnerThicknessMin = passiveInner.thickness(); + double passiveInnerThicknessMax = passiveInnerMax.thickness(); + double passiveOuterThickness = passiveOuter.thickness(); + double passiveGlueThickness = passiveGlue.thickness(); + double passiveThickness = passiveInnerThicknessMin + passiveOuterThickness + passiveGlueThickness; + double angle = passive.rotation().angle(); + + double bathRmin = caloDim.rmin(); // - margin for inclination + double bathRmax = caloDim.rmax(); // + margin for inclination + dd4hep::Tube bathOuterShape(bathRmin, bathRmax, caloDim.dz()); // make it 4 volumes + 5th for detector envelope + dd4hep::Tube bathAndServicesWarmOuterShape(cryoSideWarmDim.rmin(), cryoSideWarmDim.rmax(), caloDim.dim_z()); // make it 4 volumes + 5th for detector envelope + dd4hep::Tube bathAndServicesColdOuterShape(cryoSideColdDim.rmin(), cryoSideColdDim.rmax(), caloDim.dz()); // make it 4 volumes + 5th for detector envelope + if (cryoThicknessFront > 0) { + // 1. Create cryostat + dd4hep::Tube cryoFrontWarmShape(cryoFrontWarmDim.rmin(), cryoFrontWarmDim.rmax(), cryoFrontWarmDim.dz()); + dd4hep::Tube cryoFrontColdShape(cryoFrontColdDim.rmin(), cryoFrontColdDim.rmax(), cryoFrontColdDim.dz()); + dd4hep::Tube cryoBackColdShape(cryoBackColdDim.rmin(), cryoBackColdDim.rmax(), cryoBackColdDim.dz()); + dd4hep::Tube cryoSolenoidShape(cryoSolenoidDim.rmin(), cryoSolenoidDim.rmax(), cryoSolenoidDim.dz()); + dd4hep::Tube cryoBackWarmShape(cryoBackWarmDim.rmin(), cryoBackWarmDim.rmax(), cryoBackWarmDim.dz()); + dd4hep::Tube cryoSideColdOuterShape(cryoSideColdDim.rmin(), cryoSideColdDim.rmax(), cryoSideColdDim.dz()); + dd4hep::Tube cryoSideWarmOuterShape(cryoSideWarmDim.rmin(), cryoSideWarmDim.rmax(), cryoSideWarmDim.dz()); + dd4hep::SubtractionSolid cryoSideColdShape(cryoSideColdOuterShape, bathAndServicesColdOuterShape); + dd4hep::SubtractionSolid cryoSideWarmShape(cryoSideWarmOuterShape, bathAndServicesWarmOuterShape); + lLog << MSG::INFO << "ECAL cryostat: front warm: rmin (cm) = " << cryoFrontWarmDim.rmin() << " rmax (cm) = " << cryoFrontWarmDim.rmax() << " dz (cm) = " << cryoFrontWarmDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: front cold: rmin (cm) = " << cryoFrontColdDim.rmin() << " rmax (cm) = " << cryoFrontColdDim.rmax() << " dz (cm) = " << cryoFrontColdDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: back cold: rmin (cm) = " << cryoBackColdDim.rmin() << " rmax (cm) = " << cryoBackColdDim.rmax() << " dz (cm) = " << cryoBackColdDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: solenoid: rmin (cm) = " << cryoSolenoidDim.rmin() << " rmax (cm) = " << cryoSolenoidDim.rmax() << " dz (cm) = " << cryoSolenoidDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: back warm: rmin (cm) = " << cryoBackWarmDim.rmin() << " rmax (cm) = " << cryoBackWarmDim.rmax() << " dz (cm) = " << cryoBackWarmDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: side cold: rmin (cm) = " << cryoSideColdDim.rmin() << " rmax (cm) = " << cryoSideColdDim.rmax() << " dz (cm) = " << cryoSideColdDim.dz() - caloDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: side warm: rmin (cm) = " << cryoSideWarmDim.rmin() << " rmax (cm) = " << cryoSideWarmDim.rmax() << " dz (cm) = " << cryoSideWarmDim.dz() - caloDim.dim_z() << endmsg; + dd4hep::Volume cryoFrontWarmVol(cryostat.nameStr()+"_frontWarm", cryoFrontWarmShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoFrontColdVol(cryostat.nameStr()+"_frontCold", cryoFrontColdShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoBackColdVol(cryostat.nameStr()+"_backCold", cryoBackColdShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoSolenoidVol(cryostat.nameStr()+"_solenoid", cryoSolenoidShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoBackWarmVol(cryostat.nameStr()+"_backWarm", cryoBackWarmShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoSideWarmVol(cryostat.nameStr()+"_sideWarm", cryoSideWarmShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoSideColdVol(cryostat.nameStr()+"_sideCold", cryoSideColdShape, aLcdd.material(cryostat.materialStr())); + dd4hep::PlacedVolume cryoFrontWarmPhysVol = envelopeVol.placeVolume(cryoFrontWarmVol); + dd4hep::PlacedVolume cryoFrontColdPhysVol = envelopeVol.placeVolume(cryoFrontColdVol); + dd4hep::PlacedVolume cryoBackColdPhysVol = envelopeVol.placeVolume(cryoBackColdVol); + dd4hep::PlacedVolume cryoSolenoidPhysVol = envelopeVol.placeVolume(cryoSolenoidVol); + dd4hep::PlacedVolume cryoBackWarmPhysVol = envelopeVol.placeVolume(cryoBackWarmVol); + dd4hep::PlacedVolume cryoSideWarmPhysVol = envelopeVol.placeVolume(cryoSideWarmVol); + dd4hep::PlacedVolume cryoSideColdPhysVol = envelopeVol.placeVolume(cryoSideColdVol); + if (cryoFrontSensitive) { + cryoFrontWarmVol.setSensitiveDetector(aSensDet); + cryoFrontWarmPhysVol.addPhysVolID("cryo", 1); + cryoFrontWarmPhysVol.addPhysVolID("type", 1); + lLog << MSG::INFO << "Cryostat front warm volume set as sensitive" << endmsg; + cryoFrontColdVol.setSensitiveDetector(aSensDet); + cryoFrontColdPhysVol.addPhysVolID("cryo", 1); + cryoFrontColdPhysVol.addPhysVolID("type", 2); + lLog << MSG::INFO << "Cryostat front cold volume set as sensitive" << endmsg; + } + if (cryoBackSensitive) { + cryoBackColdVol.setSensitiveDetector(aSensDet); + cryoBackColdPhysVol.addPhysVolID("cryo", 1); + cryoBackColdPhysVol.addPhysVolID("type", 3); + lLog << MSG::INFO << "Cryostat back cold volume set as sensitive" << endmsg; + cryoSolenoidVol.setSensitiveDetector(aSensDet); + cryoSolenoidPhysVol.addPhysVolID("cryo", 1); + cryoSolenoidPhysVol.addPhysVolID("type", 4); + lLog << MSG::INFO << "Cryostat solenoid volume set as sensitive" << endmsg; + cryoBackWarmVol.setSensitiveDetector(aSensDet); + cryoBackWarmPhysVol.addPhysVolID("cryo", 1); + cryoBackWarmPhysVol.addPhysVolID("type", 5); + lLog << MSG::INFO << "Cryostat back warm volume set as sensitive" << endmsg; + } + if (cryoSideSensitive) { + cryoSideWarmVol.setSensitiveDetector(aSensDet); + cryoSideWarmPhysVol.addPhysVolID("cryo", 1); + cryoSideWarmPhysVol.addPhysVolID("type", 6); + lLog << MSG::INFO << "Cryostat side warm volume set as sensitive" << endmsg; + cryoSideColdVol.setSensitiveDetector(aSensDet); + cryoSideColdPhysVol.addPhysVolID("cryo", 1); + cryoSideColdPhysVol.addPhysVolID("type", 7); + lLog << MSG::INFO << "Cryostat side cold volume set as sensitive" << endmsg; + } + dd4hep::DetElement cryoFrontWarmDetElem(caloDetElem, "cryo_front_warm", 0); + cryoFrontWarmDetElem.setPlacement(cryoFrontWarmPhysVol); + dd4hep::DetElement cryoFrontColdDetElem(caloDetElem, "cryo_front_cold", 0); + cryoFrontColdDetElem.setPlacement(cryoFrontColdPhysVol); + dd4hep::DetElement cryoBackColdDetElem(caloDetElem, "cryo_back_cold", 0); + cryoBackColdDetElem.setPlacement(cryoBackColdPhysVol); + dd4hep::DetElement cryoSolenoidDetElem(caloDetElem, "cryo_solenoid", 0); + cryoSolenoidDetElem.setPlacement(cryoSolenoidPhysVol); + dd4hep::DetElement cryoBackWarmDetElem(caloDetElem, "cryo_back_warm", 0); + cryoBackWarmDetElem.setPlacement(cryoBackWarmPhysVol); + dd4hep::DetElement cryoSideWarmDetElem(caloDetElem, "cryo_side_warm", 0); + cryoSideWarmDetElem.setPlacement(cryoSideWarmPhysVol); + dd4hep::DetElement cryoSideColdDetElem(caloDetElem, "cryo_side_cold", 0); + cryoSideColdDetElem.setPlacement(cryoSideColdPhysVol); + // 1.2. Create place-holder for services + dd4hep::Tube servicesFrontShape(cryoFrontColdDim.rmax(), bathRmin, caloDim.dz()); + dd4hep::Tube servicesBackShape(bathRmax, cryoBackColdDim.rmin(), caloDim.dz()); + lLog << MSG::INFO << "ECAL services: front: rmin (cm) = " << cryoFrontColdDim.rmax() << " rmax (cm) = " << bathRmin << " dz (cm) = " << caloDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL services: back: rmin (cm) = " << bathRmax << " rmax (cm) = " << cryoBackColdDim.rmin() << " dz (cm) = " << caloDim.dz() << endmsg; + dd4hep::Volume servicesFrontVol("services_front", servicesFrontShape, aLcdd.material(activeMaterial)); + dd4hep::Volume servicesBackVol("services_back", servicesBackShape, aLcdd.material(activeMaterial)); + dd4hep::PlacedVolume servicesFrontPhysVol = envelopeVol.placeVolume(servicesFrontVol); + dd4hep::PlacedVolume servicesBackPhysVol = envelopeVol.placeVolume(servicesBackVol); + if (cryoFrontSensitive) { + servicesFrontVol.setSensitiveDetector(aSensDet); + servicesFrontPhysVol.addPhysVolID("cryo", 1); + servicesFrontPhysVol.addPhysVolID("type", 8); + lLog << MSG::INFO << "Services front volume set as sensitive" << endmsg; + } + if (cryoBackSensitive) { + servicesBackVol.setSensitiveDetector(aSensDet); + servicesBackPhysVol.addPhysVolID("cryo", 1); + servicesBackPhysVol.addPhysVolID("type", 9); + lLog << MSG::INFO << "Services back volume set as sensitive" << endmsg; + } + dd4hep::DetElement servicesFrontDetElem(caloDetElem, "services_front", 0); + servicesFrontDetElem.setPlacement(servicesFrontPhysVol); + dd4hep::DetElement servicesBackDetElem(caloDetElem, "services_back", 0); + servicesBackDetElem.setPlacement(servicesBackPhysVol); + } + // 2. Create bath that is inside the cryostat and surrounds the detector + // Bath is filled with active material -> but not sensitive + dd4hep::Volume bathVol(activeMaterial + "_bath", bathOuterShape, aLcdd.material(activeMaterial)); + lLog << MSG::INFO << "ECAL bath: material = " << activeMaterial << " rmin (cm) = " << bathRmin + << " rmax (cm) = " << bathRmax << " thickness in front of ECal (cm) = " << caloDim.rmin() - cryoFrontColdDim.rmax() + << " thickness behind ECal (cm) = " << cryoBackColdDim.rmin() - caloDim.rmax() << endmsg; + + // 3. Create the calorimeter by placing the passive material, trapezoid active layers, readout and again trapezoid + // active layers in the bath. + // sensitive detector for the layers + dd4hep::SensitiveDetector sd = aSensDet; + dd4hep::xml::Dimension sdType = xmlDetElem.child(_U(sensitive)); + sd.setType(sdType.typeStr()); + lLog << MSG::INFO << "ECAL calorimeter volume rmin (cm) = " << caloDim.rmin() << " rmax (cm) = " << caloDim.rmax() + << endmsg; + + // 3.a. Create the passive planes, readout in between of 2 passive planes and the remaining space fill with active + // material + ////////////////////////////// + // PASSIVE PLANES + ////////////////////////////// + lLog << MSG::INFO << "passive inner material = " << passiveInnerMaterial << "\n" + << " and outer material = " << passiveOuterMaterial << "\n" + << " thickness of inner part at inner radius (cm) = " << passiveInnerThicknessMin << "\n" + << " thickness of inner part at outer radius (cm) = " << passiveInnerThicknessMax << "\n" + << " thickness of outer part (cm) = " << passiveOuterThickness << "\n" + << " thickness of total (cm) = " << passiveThickness << "\n" + << " rotation angle = " << angle << endmsg; + uint numPlanes = + round(M_PI / asin((passiveThickness + activeThickness + readoutThickness) / (2. * caloDim.rmin() * cos(angle)))); + double dPhi = 2. * M_PI / numPlanes; + lLog << MSG::INFO << "number of passive plates = " << numPlanes << " azim. angle difference = " << dPhi << endmsg; + lLog << MSG::INFO << " distance at inner radius (cm) = " << 2. * M_PI * caloDim.rmin() / numPlanes << "\n" + << " distance at outer radius (cm) = " << 2. * M_PI * caloDim.rmax() / numPlanes << endmsg; + // Readout is in the middle between two passive planes + double offsetPassivePhi = caloDim.offset() + dPhi / 2.; + double offsetReadoutPhi = caloDim.offset() + 0; + lLog << MSG::INFO << "readout material = " << readoutMaterial << "\n" + << " thickness of readout planes (cm) = " << readoutThickness << "\n number of readout layers = " << numLayers + << endmsg; + double Rmin = caloDim.rmin(); + double Rmax = caloDim.rmax(); + double dR = Rmax - Rmin; + double planeLength = -Rmin * cos(angle) + sqrt(pow(Rmax, 2) - pow(Rmin * sin(angle), 2)); + lLog << MSG::INFO << "thickness of calorimeter (cm) = " << dR << "\n" + << " length of passive or readout planes (cm) = " << planeLength << endmsg; + + + // fill the thickness in the boundary of each layer + std::vector passiveInnerThicknessLayer(numLayers+1); + double runningHeight = 0; + for (uint iLay = 0; iLay < numLayers; iLay++) { + passiveInnerThicknessLayer[iLay] = passiveInnerThicknessMin + (passiveInnerThicknessMax - passiveInnerThicknessMin) * + (runningHeight) / (Rmax - Rmin); + runningHeight += layerHeight[iLay]; + } + passiveInnerThicknessLayer[numLayers] = passiveInnerThicknessMin + (passiveInnerThicknessMax - passiveInnerThicknessMin) * + (runningHeight) / (Rmax - Rmin); + + double passiveAngle = atan2((passiveInnerThicknessMax - passiveInnerThicknessMin) / 2., planeLength); + double cosPassiveAngle = cos(passiveAngle); + double rotatedOuterThickness = passiveOuterThickness / cosPassiveAngle; + double rotatedGlueThickness = passiveGlueThickness / cosPassiveAngle; + + // rescale layer thicknesses + double scaleLayerThickness = planeLength / layersTotalHeight; + layersTotalHeight = 0; + for (uint iLay = 0; iLay < numLayers; iLay++) { + layerHeight[iLay] *= scaleLayerThickness; + + layersTotalHeight += layerHeight[iLay]; + lLog << MSG::DEBUG << "Thickness of layer " << iLay << " : " << layerHeight[iLay] << endmsg; + } + double layerFirstOffset = -planeLength / 2. + layerHeight[0] / 2.; + + //dd4hep::Box passiveShape(passiveThickness / 2., caloDim.dz(), planeLength / 2.); + dd4hep::Trd1 passiveShape(passiveInnerThicknessMin / 2. + rotatedOuterThickness / 2. + rotatedGlueThickness / 2., + passiveInnerThicknessMax / 2. + rotatedOuterThickness / 2. + rotatedGlueThickness / 2., + caloDim.dz(), planeLength / 2.); + // inner layer is not in the first calo layer (to sample more uniformly in the layer where upstream correction is + // applied) + //dd4hep::Box passiveInnerShape(passiveInnerThickness / 2., caloDim.dz(), planeLength / 2. - layerHeight[0] / 2.); + dd4hep::Trd1 passiveInnerShape(passiveInnerThicknessLayer[1] / 2., passiveInnerThicknessMax / 2., caloDim.dz(), planeLength / 2. - layerHeight[0] / 2.); + //dd4hep::Box passiveInnerShapeFirstLayer(passiveInnerThickness / 2., caloDim.dz(), layerHeight[0] / 2.); + dd4hep::Trd1 passiveInnerShapeFirstLayer(passiveInnerThicknessMin / 2., passiveInnerThicknessLayer[1] / 2., caloDim.dz(), layerHeight[0] / 2.); + dd4hep::Box passiveOuterShape(passiveOuterThickness / 4., caloDim.dz(), planeLength / 2. / cosPassiveAngle); + dd4hep::Box passiveGlueShape(passiveGlueThickness / 4., caloDim.dz(), planeLength / 2. / cosPassiveAngle); + // passive volume consists of inner part and two outer, joind by glue + dd4hep::Volume passiveVol("passive", passiveShape, aLcdd.material("Air")); + dd4hep::Volume passiveInnerVol(passiveInnerMaterial + "_passive", passiveInnerShape, + aLcdd.material(passiveInnerMaterial)); + dd4hep::Volume passiveInnerVolFirstLayer(activeMaterial + "_passive", passiveInnerShapeFirstLayer, + aLcdd.material(activeMaterial)); + dd4hep::Volume passiveOuterVol(passiveOuterMaterial + "_passive", passiveOuterShape, + aLcdd.material(passiveOuterMaterial)); + dd4hep::Volume passiveGlueVol(passiveGlueMaterial + "_passive", passiveGlueShape, + aLcdd.material(passiveGlueMaterial)); + + if (passiveInner.isSensitive()) { + lLog << MSG::DEBUG << "Passive inner volume set as sensitive" << endmsg; + // inner part starts at second layer + double layerOffset = layerFirstOffset + layerHeight[1] / 2.; + for (uint iLayer = 1; iLayer < numLayers; iLayer++) { + //dd4hep::Box layerPassiveInnerShape(passiveInnerThickness / 2., caloDim.dz(), layerHeight[iLayer] / 2.); + dd4hep::Trd1 layerPassiveInnerShape(passiveInnerThicknessLayer[iLayer] / 2., passiveInnerThicknessLayer[iLayer+1] / 2., caloDim.dz(), layerHeight[iLayer] / 2.); + dd4hep::Volume layerPassiveInnerVol(passiveInnerMaterial, layerPassiveInnerShape, + aLcdd.material(passiveInnerMaterial)); + layerPassiveInnerVol.setSensitiveDetector(aSensDet); + dd4hep::PlacedVolume layerPassiveInnerPhysVol = + passiveInnerVol.placeVolume(layerPassiveInnerVol, dd4hep::Position(0, 0, layerOffset)); + layerPassiveInnerPhysVol.addPhysVolID("layer", iLayer); + dd4hep::DetElement layerPassiveInnerDetElem("layer", iLayer); + layerPassiveInnerDetElem.setPlacement(layerPassiveInnerPhysVol); + if (iLayer != numLayers - 1) { + layerOffset += layerHeight[iLayer] / 2. + layerHeight[iLayer + 1] / 2.; + } + } + } + if (passiveOuter.isSensitive()) { + lLog << MSG::DEBUG << "Passive outer volume set as sensitive" << endmsg; + double layerOffset = layerFirstOffset / cosPassiveAngle; + for (uint iLayer = 0; iLayer < numLayers; iLayer++) { + dd4hep::Box layerPassiveOuterShape(passiveOuterThickness / 4., caloDim.dz(), layerHeight[iLayer] / 2. / cosPassiveAngle); + dd4hep::Volume layerPassiveOuterVol(passiveOuterMaterial, layerPassiveOuterShape, + aLcdd.material(passiveOuterMaterial)); + layerPassiveOuterVol.setSensitiveDetector(aSensDet); + dd4hep::PlacedVolume layerPassiveOuterPhysVol = + passiveOuterVol.placeVolume(layerPassiveOuterVol, dd4hep::Position(0, 0, layerOffset)); + layerPassiveOuterPhysVol.addPhysVolID("layer", iLayer); + dd4hep::DetElement layerPassiveOuterDetElem("layer", iLayer); + layerPassiveOuterDetElem.setPlacement(layerPassiveOuterPhysVol); + if (iLayer != numLayers - 1) { + layerOffset += (layerHeight[iLayer] / 2. + layerHeight[iLayer + 1] / 2.) / cosPassiveAngle; + } + } + } + if (passiveGlue.isSensitive()) { + lLog << MSG::DEBUG << "Passive glue volume set as sensitive" << endmsg; + double layerOffset = layerFirstOffset / cosPassiveAngle; + for (uint iLayer = 0; iLayer < numLayers; iLayer++) { + dd4hep::Box layerPassiveGlueShape(passiveGlueThickness / 4., caloDim.dz(), layerHeight[iLayer] / 2. / cosPassiveAngle); + dd4hep::Volume layerPassiveGlueVol(passiveGlueMaterial, layerPassiveGlueShape, + aLcdd.material(passiveGlueMaterial)); + layerPassiveGlueVol.setSensitiveDetector(aSensDet); + dd4hep::PlacedVolume layerPassiveGluePhysVol = + passiveGlueVol.placeVolume(layerPassiveGlueVol, dd4hep::Position(0, 0, layerOffset)); + layerPassiveGluePhysVol.addPhysVolID("layer", iLayer); + dd4hep::DetElement layerPassiveGlueDetElem("layer", iLayer); + layerPassiveGlueDetElem.setPlacement(layerPassiveGluePhysVol); + if (iLayer != numLayers - 1) { + layerOffset += (layerHeight[iLayer] / 2. + layerHeight[iLayer + 1] / 2.) / cosPassiveAngle; + } + } + } + + dd4hep::PlacedVolume passiveInnerPhysVol = + passiveVol.placeVolume(passiveInnerVol, dd4hep::Position(0, 0, layerHeight[0] / 2.)); + dd4hep::PlacedVolume passiveInnerPhysVolFirstLayer = + passiveVol.placeVolume(passiveInnerVolFirstLayer, dd4hep::Position(0, 0, layerFirstOffset)); + dd4hep::PlacedVolume passiveOuterPhysVolBelow = passiveVol.placeVolume( + passiveOuterVol, + dd4hep::Transform3D(dd4hep::RotationY(-passiveAngle), + dd4hep::Position(-(passiveInnerThicknessMin + passiveInnerThicknessMax) / 4. - + rotatedGlueThickness / 2. - rotatedOuterThickness / 4., 0, 0))); + dd4hep::PlacedVolume passiveOuterPhysVolAbove = passiveVol.placeVolume( + passiveOuterVol, + dd4hep::Transform3D(dd4hep::RotationY(passiveAngle), + dd4hep::Position((passiveInnerThicknessMin + passiveInnerThicknessMax) / 4. + + rotatedGlueThickness / 2. + rotatedOuterThickness / 4., 0, 0))); + dd4hep::PlacedVolume passiveGluePhysVolBelow = passiveVol.placeVolume( + passiveGlueVol, + dd4hep::Transform3D(dd4hep::RotationY(-passiveAngle), + dd4hep::Position(-(passiveInnerThicknessMin + passiveInnerThicknessMax) / 4. - + rotatedGlueThickness / 4., 0, 0))); + dd4hep::PlacedVolume passiveGluePhysVolAbove = passiveVol.placeVolume( + passiveGlueVol, + dd4hep::Transform3D(dd4hep::RotationY(passiveAngle), + dd4hep::Position((passiveInnerThicknessMin + passiveInnerThicknessMax) / 4. + + rotatedGlueThickness / 4., 0, 0))); + passiveInnerPhysVol.addPhysVolID("subtype", 0); + passiveInnerPhysVolFirstLayer.addPhysVolID("subtype", 0); + passiveOuterPhysVolBelow.addPhysVolID("subtype", 1); + passiveOuterPhysVolAbove.addPhysVolID("subtype", 2); + passiveGluePhysVolBelow.addPhysVolID("subtype", 3); + passiveGluePhysVolAbove.addPhysVolID("subtype", 4); + if (passiveInner.isSensitive()) { + passiveInnerVolFirstLayer.setSensitiveDetector(aSensDet); + passiveInnerPhysVolFirstLayer.addPhysVolID("layer", 0); + dd4hep::DetElement passiveInnerDetElemFirstLayer("layer", 0); + passiveInnerDetElemFirstLayer.setPlacement(passiveInnerPhysVolFirstLayer); + } + + ////////////////////////////// + // READOUT PLANES + ////////////////////////////// + dd4hep::Box readoutShape(readoutThickness / 2., caloDim.dz(), planeLength / 2.); + dd4hep::Volume readoutVol(readoutMaterial, readoutShape, aLcdd.material(readoutMaterial)); + if (readout.isSensitive()) { + lLog << MSG::INFO << "Readout volume set as sensitive" << endmsg; + double layerOffset = layerFirstOffset; + for (uint iLayer = 0; iLayer < numLayers; iLayer++) { + dd4hep::Box layerReadoutShape(readoutThickness / 2., caloDim.dz(), layerHeight[iLayer] / 2.); + dd4hep::Volume layerReadoutVol(readoutMaterial, layerReadoutShape, aLcdd.material(readoutMaterial)); + layerReadoutVol.setSensitiveDetector(aSensDet); + dd4hep::PlacedVolume layerReadoutPhysVol = + readoutVol.placeVolume(layerReadoutVol, dd4hep::Position(0, 0, layerOffset)); + layerReadoutPhysVol.addPhysVolID("layer", iLayer); + dd4hep::DetElement layerReadoutDetElem("layer", iLayer); + layerReadoutDetElem.setPlacement(layerReadoutPhysVol); + if (iLayer != numLayers - 1) { + layerOffset += layerHeight[iLayer] / 2. + layerHeight[iLayer + 1] / 2.; + } + } + } + + ////////////////////////////// + // ACTIVE + ////////////////////////////// + // thickness of active layers at inner radius and outer ( = distance between passive plane and readout plane) + // at inner radius: distance projected at plane perpendicular to readout plane + double activeInThickness = Rmin * sin(dPhi / 2.) * cos(angle); + activeInThickness -= passiveThickness * (0.5 - activePassiveOverlap); + // at outer radius: distance projected at plane perpendicular to readout plane + double activeOutThickness = (Rmin + planeLength) * sin(dPhi / 2.) * cos(angle); + // make correction for outer readius caused by inclination angle + // first calculate intersection of readout plane and plane parallel to shifted passive plane + double xIntersect = (Rmin * (tan(angle) - cos(dPhi / 2.) * tan(angle + dPhi / 2.)) - planeLength * sin(dPhi / 2.)) / + (tan(angle) - tan(angle + dPhi / 2.)); + double yIntersect = tan(angle) * xIntersect + Rmin * (sin(dPhi / 2.) - tan(angle)) + planeLength * sin(dPhi / 2.); + // distance from inner radius to intersection + double correction = + planeLength - sqrt(pow(xIntersect - Rmin * cos(dPhi / 2), 2) + pow(yIntersect - Rmin * sin(dPhi / 2), 2)); + // correction to the active thickness + activeOutThickness += 2. * correction * sin(dPhi / 4.); + activeOutThickness -= passiveThickness * (0.5 - activePassiveOverlap); + // print the active layer dimensions + double activeInThicknessAfterSubtraction = + 2. * activeInThickness - readoutThickness - 2. * activePassiveOverlap * passiveThickness; + double activeOutThicknessAfterSubtraction = + 2. * activeOutThickness - readoutThickness - 2. * activePassiveOverlap * + (passiveThickness + passiveInnerThicknessMax - passiveInnerThicknessMin); // correct thickness for trapezoid + lLog << MSG::INFO << "active material = " << activeMaterial + << " active layers thickness at inner radius (cm) = " << activeInThicknessAfterSubtraction + << " thickness at outer radious (cm) = " << activeOutThicknessAfterSubtraction << " making " + << (activeOutThicknessAfterSubtraction - activeInThicknessAfterSubtraction) * 100 / + activeInThicknessAfterSubtraction + << " % increase." << endmsg; + lLog << MSG::INFO + << "active passive initial overlap (before subtraction) (cm) = " << passiveThickness * activePassiveOverlap + << " = " << activePassiveOverlap * 100 << " %" << endmsg; + + // creating shape for rows of layers (active material between two passive planes, with readout in the middle) + // first define area between two passive planes, area can reach up to the symmetry axis of passive plane + dd4hep::Trd1 activeOuterShape(activeInThickness, activeOutThickness, caloDim.dz(), planeLength / 2.); + // subtract readout shape from the middle + dd4hep::SubtractionSolid activeShapeNoReadout(activeOuterShape, readoutShape); + + // make calculation for active plane that is inclined with 0 deg (= offset + angle) + double Cx = Rmin * cos(-angle) + planeLength / 2.; + double Cy = Rmin * sin(-angle); + double Ax = Rmin * cos(-angle + dPhi / 2.) + planeLength / 2. * cos(dPhi / 2.); + double Ay = Rmin * sin(-angle + dPhi / 2.) + planeLength / 2. * sin(dPhi / 2.); + double CAx = fabs(Ax - Cx); + double CAy = fabs(Ay - Cy); + double zprim, xprim; + zprim = CAx; + xprim = CAy; + + double Bx = Rmin * cos(-angle - dPhi / 2.) + planeLength / 2. * cos(-dPhi / 2.); + double By = Rmin * sin(-angle - dPhi / 2.) + planeLength / 2. * sin(-dPhi / 2.); + double CBx = fabs(Bx - Cx); + double CBy = fabs(By - Cy); + double zprimB, xprimB; + zprimB = CBx; + xprimB = CBy; + + // subtract passive volume above + dd4hep::SubtractionSolid activeShapeNoPassiveAbove( + activeShapeNoReadout, passiveShape, + dd4hep::Transform3D(dd4hep::RotationY(-dPhi / 2.), + dd4hep::Position(-fabs(xprim), 0, fabs(zprim)))); + // subtract passive volume below + dd4hep::SubtractionSolid activeShape( + activeShapeNoPassiveAbove, passiveShape, + dd4hep::Transform3D(dd4hep::RotationY(dPhi / 2.), + dd4hep::Position(fabs(xprimB), 0, -fabs(zprimB)))); + dd4hep::Volume activeVol("active", activeShape, aLcdd.material("Air")); + + std::vector layerPhysVols; + // place layers within active volume + std::vector layerInThickness; + std::vector layerOutThickness; + double layerIncreasePerUnitThickness = (activeOutThickness - activeInThickness) / layersTotalHeight; + for (uint iLay = 0; iLay < numLayers; iLay++) { + if (iLay == 0) { + layerInThickness.push_back(activeInThickness); + } else { + layerInThickness.push_back(layerOutThickness[iLay - 1]); + } + layerOutThickness.push_back(layerInThickness[iLay] + layerIncreasePerUnitThickness * layerHeight[iLay]); + } + double layerOffset = layerFirstOffset; + for (uint iLayer = 0; iLayer < numLayers; iLayer++) { + dd4hep::Trd1 layerOuterShape(layerInThickness[iLayer], layerOutThickness[iLayer], caloDim.dz(), layerHeight[iLayer] / 2.); + dd4hep::SubtractionSolid layerShapeNoReadout(layerOuterShape, readoutShape); + dd4hep::SubtractionSolid layerShapeNoPassiveAbove( + layerShapeNoReadout, passiveShape, + dd4hep::Transform3D(dd4hep::RotationY(-dPhi / 2.), + dd4hep::Position(-fabs(xprim), 0, fabs(zprim) - layerOffset))); + // subtract passive volume below + dd4hep::SubtractionSolid layerShape( + layerShapeNoPassiveAbove, passiveShape, + dd4hep::Transform3D(dd4hep::RotationY(dPhi / 2.), + dd4hep::Position(fabs(xprimB), 0, -fabs(zprimB) - layerOffset))); + dd4hep::Volume layerVol("layer", layerShape, aLcdd.material(activeMaterial)); + layerVol.setSensitiveDetector(aSensDet); + layerPhysVols.push_back(activeVol.placeVolume(layerVol, dd4hep::Position(0, 0, layerOffset))); + layerPhysVols.back().addPhysVolID("layer", iLayer); + if (iLayer != numLayers - 1) { + layerOffset += layerHeight[iLayer] / 2. + layerHeight[iLayer + 1] / 2.; + } + } + + dd4hep::DetElement bathDetElem(caloDetElem, "bath", 1); + std::vector activePhysVols; + // Next place elements: passive planes, readout planes and rows of layers + for (uint iPlane = 0; iPlane < numPlanes; iPlane++) { + // first calculate positions of passive and readout planes + // PASSIVE + // calculate centre position of the plane without plane rotation + double phi = offsetPassivePhi + iPlane * dPhi; + double xRadial = (Rmin + planeLength / 2.) * cos(phi); + double yRadial = (Rmin + planeLength / 2.) * sin(phi); + // calculate position of the beginning of plane + double xRmin = Rmin * cos(phi); + double yRmin = Rmin * sin(phi); + // rotate centre by angle wrt beginning of plane + double xRotated = xRmin + (xRadial - xRmin) * cos(angle) - (yRadial - yRmin) * sin(angle); + double yRotated = yRmin + (xRadial - xRmin) * sin(angle) + (yRadial - yRmin) * cos(angle); + dd4hep::Transform3D transform(dd4hep::RotationX(-M_PI / 2.) // to get in XY plane + * + dd4hep::RotationY(M_PI / 2. // to get pointed towards centre + - + phi - angle), + dd4hep::Position(xRotated, yRotated, 0)); + dd4hep::PlacedVolume passivePhysVol = bathVol.placeVolume(passiveVol, transform); + passivePhysVol.addPhysVolID("module", iPlane); + passivePhysVol.addPhysVolID("type", 1); // 0 = active, 1 = passive, 2 = readout + dd4hep::DetElement passiveDetElem(bathDetElem, "passive" + std::to_string(iPlane), iPlane); + passiveDetElem.setPlacement(passivePhysVol); + + // READOUT + // calculate centre position of the plane without plane rotation + double phiRead = offsetReadoutPhi + iPlane * dPhi; + double xRadialRead = (Rmin + planeLength / 2.) * cos(phiRead); + double yRadialRead = (Rmin + planeLength / 2.) * sin(phiRead); + // calculate position of the beginning of plane + double xRminRead = Rmin * cos(phiRead); + double yRminRead = Rmin * sin(phiRead); + // rotate centre by angle wrt beginning of plane + double xRotatedRead = xRminRead + (xRadialRead - xRminRead) * cos(angle) - (yRadialRead - yRminRead) * sin(angle); + double yRotatedRead = yRminRead + (xRadialRead - xRminRead) * sin(angle) + (yRadialRead - yRminRead) * cos(angle); + dd4hep::Transform3D transformRead( + dd4hep::RotationX(-M_PI / 2.) // to get in XY plane + * + dd4hep::RotationY(M_PI / 2. // to get pointed towards centre + - + phiRead - angle), + dd4hep::Position(xRotatedRead, yRotatedRead, 0)); + dd4hep::PlacedVolume readoutPhysVol = bathVol.placeVolume(readoutVol, transformRead); + readoutPhysVol.addPhysVolID("module", iPlane); + readoutPhysVol.addPhysVolID("type", 2); // 0 = active, 1 = passive, 2 = readout + dd4hep::DetElement readoutDetElem(bathDetElem, "readout" + std::to_string(iPlane), iPlane); + readoutDetElem.setPlacement(readoutPhysVol); + + // ACTIVE + dd4hep::Rotation3D rotationActive(dd4hep::RotationX(-M_PI / 2) * + dd4hep::RotationY(M_PI / 2 - phiRead - angle)); + activePhysVols.push_back(bathVol.placeVolume( + activeVol, + dd4hep::Transform3D(rotationActive, dd4hep::Position(xRotatedRead, yRotatedRead, 0)))); + activePhysVols.back().addPhysVolID("module", iPlane); + activePhysVols.back().addPhysVolID("type", 0); // 0 = active, 1 = passive, 2 = readout + } + dd4hep::PlacedVolume bathPhysVol = envelopeVol.placeVolume(bathVol); + bathDetElem.setPlacement(bathPhysVol); + for (uint iPlane = 0; iPlane < numPlanes; iPlane++) { + dd4hep::DetElement activeDetElem(bathDetElem, "active" + std::to_string(iPlane), iPlane); + activeDetElem.setPlacement(activePhysVols[iPlane]); + for (uint iLayer = 0; iLayer < numLayers; iLayer++) { + dd4hep::DetElement layerDetElem(activeDetElem, "layer" + std::to_string(iLayer), iLayer); + layerDetElem.setPlacement(layerPhysVols[iLayer]); + } + } + + // Place the envelope + dd4hep::Volume motherVol = aLcdd.pickMotherVolume(caloDetElem); + dd4hep::PlacedVolume envelopePhysVol = motherVol.placeVolume(envelopeVol); + envelopePhysVol.addPhysVolID("system", xmlDetElem.id()); + caloDetElem.setPlacement(envelopePhysVol); + + // Create caloData object + auto caloData = new dd4hep::rec::LayeredCalorimeterData; + caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::BarrelLayout; + caloDetElem.addExtension(caloData); + + // Set type flags + dd4hep::xml::setDetectorTypeFlag(xmlDetElem, caloDetElem); + + return caloDetElem; +} +} // namespace det + +DECLARE_DETELEMENT(EmCaloBarrelInclinedFCCee, det::createECalBarrelInclinedFCCee)