diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml new file mode 100644 index 00000000..49c76c53 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml @@ -0,0 +1,46 @@ + + + + + + Master compact file describing the latest developments of the FCCee IDEA detector concept with a LAr calorimeter. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/BeamInstrumentation.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/BeamInstrumentation.xml new file mode 100644 index 00000000..49647b2b --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/BeamInstrumentation.xml @@ -0,0 +1,36 @@ + + + + COmpensating and screening solenoids for FCCee + + + + + Beampipe Instrumentation + + + + + + +
+ + + + + + + + +
+ +
+ + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Beampipe.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Beampipe.xml new file mode 100644 index 00000000..b86af521 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Beampipe.xml @@ -0,0 +1,150 @@ + + + + A beampipe for FCCee, R(central) = 1.5 cm + + + + + + + + + + + + + + + + + + + + Part of beampipe made of Beryllium + + + + + + + +
+ +
+ + + + + + + + + + + + + Golden foil in the inner part of the Be beampipe + +
+ +
+ + Part of beampipe made of Copper + +
+ + + + + +
+ + + +
+ + +
+ + +Full Cone Tungsten Shield + + + + Before HOM space +
+ + After HOM space (1197.5*m - 1298.7*mm) +18 cm as solenoid is now closer to IP +
+ + +Asymmetric Tungsten Shield no Rotation + + + + +
+ + was 370 +
+ + one degree less, to fit lumical window +
+ +
+ + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_DectDimensions.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_DectDimensions.xml new file mode 100644 index 00000000..439e20cd --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_DectDimensions.xml @@ -0,0 +1,230 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_DectEmptyMaster.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_DectEmptyMaster.xml new file mode 100644 index 00000000..a5ad3e23 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_DectEmptyMaster.xml @@ -0,0 +1,30 @@ + + + + + + + + + + + + Use this one if you want to use official dimensions but only place one detector inside + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged.xml new file mode 100644 index 00000000..4fd78287 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged.xml @@ -0,0 +1,144 @@ + + + + + + Settings for the inclined EM calorimeter. + 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:3,subtype:3,layer:8,module:11,theta:10 + + + + + + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_EcalEndcaps_coneCryo.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_EcalEndcaps_coneCryo.xml new file mode 100644 index 00000000..33abcf03 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_EcalEndcaps_coneCryo.xml @@ -0,0 +1,167 @@ + + + + + + Liquid argon EM calorimeter endcap design. + Electromagnetic part (EMEC) includes lead+steel absorber. + Absorbers have shape of simple discs. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,subsystem:1,type:3,subtype:3,layer:8,sublayer:8,eta:10,phi:10 + + + + system:4,subsystem:1,type:3,subtype:3,layer:8,eta:10,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml new file mode 100644 index 00000000..0acecdb1 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml @@ -0,0 +1,133 @@ + + + + The first FCCee HCal layout based on ATLAS HCal, not optimised yet + 1. Nov 2022, J. Faltova: update material and radial segmentation for FCCee + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,layer:5,row:9,eta:9,phi:10 + + + + system:4,layer:5,eta:9,phi:10 + + + + + system:4,layer:5,row:9,eta:0,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalEndcaps_ThreeParts_TileCal.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalEndcaps_ThreeParts_TileCal.xml new file mode 100644 index 00000000..c5cbe72a --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalEndcaps_ThreeParts_TileCal.xml @@ -0,0 +1,141 @@ + + + + HCal layout based on ATLAS HCal, with realistic longitudinal segmentation and steel support + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,type:2,layer:5,row:9,eta:11,phi:10 + + + + system:4,type:2,layer:4,eta:11,phi:10 + + + + + system:4,type:2,layer:5,row:9,eta:10,phi:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/HOMAbsorber.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/HOMAbsorber.xml new file mode 100644 index 00000000..862207ea --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/HOMAbsorber.xml @@ -0,0 +1,64 @@ + + + + Higher mode absorber for FCCee + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/LumiCal.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/LumiCal.xml new file mode 100644 index 00000000..15f0434e --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/LumiCal.xml @@ -0,0 +1,191 @@ + + + + LumiCal for FCCee detector based on CLD + + + + + + + + + + system:8,barrel:3,layer:8,slice:8,r:32:-16,phi:-16 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/MuonTagger.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/MuonTagger.xml new file mode 100644 index 00000000..79b0ed18 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/MuonTagger.xml @@ -0,0 +1,50 @@ + + + + + + Simple muon tagger - barrel and endcaps + + + + + + + + + + + system:4,subsystem:1,type:3,subtype:3,layer:8,sublayer:8,eta:10,phi:10 + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/SimplifiedDriftChamber.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/SimplifiedDriftChamber.xml new file mode 100644 index 00000000..11586a5d --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/SimplifiedDriftChamber.xml @@ -0,0 +1,58 @@ + + + + + + A simplified implementation of the drift chamber for the FCCee-IDEA concept + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ${GlobalTrackerReadoutID_DCH} + + + + + + + + Dimensions for the drift chamber + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Solenoid_o1_v01_02.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Solenoid_o1_v01_02.xml new file mode 100644 index 00000000..69a117f8 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Solenoid_o1_v01_02.xml @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + Solenoid + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Vertex.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Vertex.xml new file mode 100644 index 00000000..a28449d6 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/Vertex.xml @@ -0,0 +1,234 @@ + + + + + + CLD Vertex Detector for FCCee + + + Tracking detectors + + + + + + + + + + + + Vertex Assembly + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ${GlobalTrackerReadoutID} + + + ${GlobalTrackerReadoutID} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Vertex Detector Endcaps + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/elements.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/elements.xml new file mode 100644 index 00000000..f35eb345 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/elements.xml @@ -0,0 +1,884 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/materials.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/materials.xml new file mode 100644 index 00000000..a899266c --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/materials.xml @@ -0,0 +1,261 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FCCee/ALLEGRO/compact/README.md b/FCCee/ALLEGRO/compact/README.md index ba8b5c28..9a6ded62 100644 --- a/FCCee/ALLEGRO/compact/README.md +++ b/FCCee/ALLEGRO/compact/README.md @@ -1,3 +1,5 @@ ALLEGRO ======================== -ALLEGRO_o1_v01: it is a liquid Noble gas based detector. This version picked from the latest version in FCCDetectors repo. +ALLEGRO_o1_v01: it is a liquid Noble gas based detector. This version picked from the latest version in FCCDetectors repo. + +ALLEGRO_o1_v02: evolves from o1_v01, replacing the barrel ECAL. This version has a constant cell size in theta for the ECAL barrel (instead of eta as in o1_v01) and now it is possible to have a different number of cells merged for each longitudinal layer. diff --git a/detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v02_geo.cpp b/detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v02_geo.cpp new file mode 100644 index 00000000..9f62dec8 --- /dev/null +++ b/detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v02_geo.cpp @@ -0,0 +1,618 @@ +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Handle.h" +#include "XML/Utilities.h" + +#include + +// Taken from https://github.com/HEP-FCC/FCCDetectors/blob/main/Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp (go there to see older commit history) + +// 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 createECalBarrelInclined(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::Dimension cryoDim(cryostat.dimensions()); + double cryoThicknessFront = cryoDim.rmin2() - cryoDim.rmin1(); + dd4hep::xml::DetElement cryoFront = cryostat.child(_Unicode(front)); + dd4hep::xml::DetElement cryoBack = cryostat.child(_Unicode(back)); + dd4hep::xml::DetElement cryoSide = cryostat.child(_Unicode(side)); + bool cryoFrontSensitive = cryoFront.isSensitive(); + bool cryoBackSensitive = cryoBack.isSensitive(); + bool cryoSideSensitive = cryoSide.isSensitive(); + + // 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 << "Number of layers: " << numLayers << " total thickness " << layersTotalHeight << endmsg; + // The following code checks if the xml geometry file contains a constant defining + // the number of layers the barrel. In that case, it makes the program abort + // if the number of planes in the xml is different from the one calculated from + // the geometry. This is because the number of layers is needed + // in other parts of the code (the readout for the FCC-ee ECAL with + // inclined modules). + int nLayers = -1; + try { + nLayers = aLcdd.constant("ECalBarrelNumLayers"); + } + catch(...) { + ; + } + if (nLayers > 0 && nLayers != int(numLayers)) { + lLog << MSG::ERROR << "Incorrect number of layers (ECalBarrelNumLayers) in xml file!" << endmsg; + // todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure")); + // make the code crash (incidentSvc does not work) + // Andre, Alvaro, assert replaced by exception + throw std::runtime_error("Incorrect number of layers (ECalBarrelNumLayers) in xml file!"); + } + + 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 bathAndServicesOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), caloDim.dz()); // make it 4 volumes + 5th for detector envelope + if (cryoThicknessFront > 0) { + // 1. Create cryostat + dd4hep::Tube cryoFrontShape(cryoDim.rmin1(), cryoDim.rmin2(), cryoDim.dz()); + dd4hep::Tube cryoBackShape(cryoDim.rmax1(), cryoDim.rmax2(), cryoDim.dz()); + dd4hep::Tube cryoSideOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz()); + dd4hep::SubtractionSolid cryoSideShape(cryoSideOuterShape, bathAndServicesOuterShape); + lLog << MSG::INFO << "ECAL cryostat: front: rmin (cm) = " << cryoDim.rmin1() << " rmax (cm) = " << cryoDim.rmin2() << " dz (cm) = " << cryoDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: back: rmin (cm) = " << cryoDim.rmax1() << " rmax (cm) = " << cryoDim.rmax2() << " dz (cm) = " << cryoDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL cryostat: side: rmin (cm) = " << cryoDim.rmin2() << " rmax (cm) = " << cryoDim.rmax1() << " dz (cm) = " << cryoDim.dz() - caloDim.dz() << endmsg; + dd4hep::Volume cryoFrontVol(cryostat.nameStr()+"_front", cryoFrontShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoBackVol(cryostat.nameStr()+"_back", cryoBackShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoSideVol(cryostat.nameStr()+"_side", cryoSideShape, aLcdd.material(cryostat.materialStr())); + dd4hep::PlacedVolume cryoFrontPhysVol = envelopeVol.placeVolume(cryoFrontVol); + dd4hep::PlacedVolume cryoBackPhysVol = envelopeVol.placeVolume(cryoBackVol); + dd4hep::PlacedVolume cryoSidePhysVol = envelopeVol.placeVolume(cryoSideVol); + if (cryoFrontSensitive) { + cryoFrontVol.setSensitiveDetector(aSensDet); + cryoFrontPhysVol.addPhysVolID("cryo", 1); + cryoFrontPhysVol.addPhysVolID("type", 1); + lLog << MSG::INFO << "Cryostat front volume set as sensitive" << endmsg; + } + if (cryoBackSensitive) { + cryoBackVol.setSensitiveDetector(aSensDet); + cryoBackPhysVol.addPhysVolID("cryo", 1); + cryoBackPhysVol.addPhysVolID("type", 2); + lLog << MSG::INFO << "Cryostat back volume set as sensitive" << endmsg; + } + if (cryoSideSensitive) { + cryoSideVol.setSensitiveDetector(aSensDet); + cryoSidePhysVol.addPhysVolID("cryo", 1); + cryoSidePhysVol.addPhysVolID("type", 3); + lLog << MSG::INFO << "Cryostat front volume set as sensitive" << endmsg; + } + dd4hep::DetElement cryoFrontDetElem(caloDetElem, "cryo_front", 0); + cryoFrontDetElem.setPlacement(cryoFrontPhysVol); + dd4hep::DetElement cryoBackDetElem(caloDetElem, "cryo_back", 0); + cryoBackDetElem.setPlacement(cryoBackPhysVol); + dd4hep::DetElement cryoSideDetElem(caloDetElem, "cryo_side", 0); + cryoSideDetElem.setPlacement(cryoSidePhysVol); + // 1.2. Create place-holder for services + dd4hep::Tube servicesFrontShape(cryoDim.rmin2(), bathRmin, caloDim.dz()); + dd4hep::Tube servicesBackShape(bathRmax, cryoDim.rmax1(), caloDim.dz()); + lLog << MSG::INFO << "ECAL services: front: rmin (cm) = " << cryoDim.rmin2() << " rmax (cm) = " << bathRmin << " dz (cm) = " << caloDim.dz() << endmsg; + lLog << MSG::INFO << "ECAL services: back: rmin (cm) = " << bathRmax << " rmax (cm) = " << cryoDim.rmax1() << " 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", 4); + lLog << MSG::INFO << "Services front volume set as sensitive" << endmsg; + } + if (cryoBackSensitive) { + servicesBackVol.setSensitiveDetector(aSensDet); + servicesBackPhysVol.addPhysVolID("cryo", 1); + servicesBackPhysVol.addPhysVolID("type", 5); + 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() - cryoDim.rmin2() + << " thickness behind ECal (cm) = " << cryoDim.rmax1() - 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; + + // The following code checks if the xml geometry file contains a constant defining + // the number of planes in the barrel. In that case, it makes the program abort + // if the number of planes in the xml is different from the one calculated from + // the geometry. This is because the number of plane information (retrieved from the + // xml) is used in other parts of the code (the readout for the FCC-ee ECAL with + // inclined modules). In principle the code above should be refactored so that the number + // of planes is one of the inputs of the calculation and other geometrical parameters + // are adjusted accordingly. This is left for the future, and we use the workaround + // below to enforce for the time being that the number of planes is "correct" + int nModules = -1; + try { + nModules = aLcdd.constant("ECalBarrelNumPlanes"); + } + catch(...) { + ; + } + if (nModules > 0 && nModules != int(numPlanes)) { + lLog << MSG::ERROR << "Incorrect number of planes (ECalBarrelNumPlanes) in xml file!" << endmsg; + // todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure")); + // make the code crash (incidentSvc does not work) + // Andre, Alvaro, assert replaced by exception + throw std::runtime_error("Incorrect number of planes (ECalBarrelNumPlanes) in xml file!"); + } + // 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(ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v02, det::createECalBarrelInclined) diff --git a/detector/calorimeter/README.md b/detector/calorimeter/README.md index d4596f46..e99fa66f 100644 --- a/detector/calorimeter/README.md +++ b/detector/calorimeter/README.md @@ -11,6 +11,9 @@ The documentation about its usage is [here](../../doc/detector/calorimeter/ECalB ### o1_v01 Original version taken from [FCCDetectors](https://github.com/HEP-FCC/FCCDetectors/blob/main/Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp). +### o1_v02 +New version adapted to the theta segmentation with the possibility to have different number of cell merged per layer. The main difference is that now one has to set `ECalBarrelNumLayers` and `ECalBarrelNumPlanes` in the xml while before it was just dynamically computed based on other xml parameters (also, to avoid silent mistakes, the number from the xml and the one computed dynamically must match). + ## CaloDisks This sub-detector makes calorimeter endcaps (original and reflected). It is used in ALLEGRO detector concept. diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h new file mode 100644 index 00000000..c0b79ea8 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h @@ -0,0 +1,127 @@ +#ifndef DD4HEP_DDCORE_GRIDMODULETHETAMERGED_H +#define DD4HEP_DDCORE_GRIDMODULETHETAMERGED_H 1 + +// FCCSW +#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" + +// DD4hep +#include "DD4hep/Segmentations.h" +#include "DD4hep/detail/SegmentationsInterna.h" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + +/// Namespace for base segmentations + + +// Forward declarations +class Segmentation; +template +class SegmentationWrapper; + +/// We need some abbreviation to make the code more readable. +typedef Handle> FCCSWGridModuleThetaMergedHandle_k4geo; + +/// Implementation class for the module-theta (with per-layer merging) segmentation. +/** + * Concrete user handle to serve specific needs of client code + * which requires access to the base functionality not served + * by the super-class Segmentation. + * + * Note: + * We only check the validity of the underlying handle. + * If for whatever reason the implementation object is not valid + * This is not checked. + * In principle this CANNOT happen unless some brain-dead has + * fiddled with the handled object directly..... + * + * Note: + * The handle base corrsponding to this object in for + * conveniance reasons instantiated in DD4hep/src/Segmentations.cpp. + * + */ +class FCCSWGridModuleThetaMerged_k4geo : public FCCSWGridModuleThetaMergedHandle_k4geo { +public: + /// Definition of the basic handled object + typedef FCCSWGridModuleThetaMergedHandle_k4geo::Object Object; + +public: + /// Default constructor + FCCSWGridModuleThetaMerged_k4geo() = default; + /// Copy constructor + FCCSWGridModuleThetaMerged_k4geo(const FCCSWGridModuleThetaMerged_k4geo& e) = default; + /// Copy Constructor from segmentation base object + FCCSWGridModuleThetaMerged_k4geo(const Segmentation& e) : Handle(e) {} + /// Copy constructor from handle + FCCSWGridModuleThetaMerged_k4geo(const Handle& e) : Handle(e) {} + /// Copy constructor from other polymorph/equivalent handle + template + FCCSWGridModuleThetaMerged_k4geo(const Handle& e) : Handle(e) {} + /// Assignment operator + FCCSWGridModuleThetaMerged_k4geo& operator=(const FCCSWGridModuleThetaMerged_k4geo& seg) = default; + /// Equality operator + bool operator==(const FCCSWGridModuleThetaMerged_k4geo& seg) const { return m_element == seg.m_element; } + + /// determine the position based on the cell ID + inline Position position(const CellID& id) const { return Position(access()->implementation->position(id)); } + + /// 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); + } + + /// access the grid size in theta + inline double gridSizeTheta() const { return access()->implementation->gridSizeTheta(); } + + /// access the coordinate offset in theta + inline double offsetTheta() const { return access()->implementation->offsetTheta(); } + + /// access the number of theta cells merged for given layer + inline int mergedThetaCells(const int layer) const { return access()->implementation->mergedThetaCells(layer); } + + /// access the number of modules + inline int nModules() const { return access()->implementation->nModules(); } + + /// access the number of layers + inline int nLayers() const { return access()->implementation->nLayers(); } + + /// access the number of merged modules for given layer + inline int mergedModules(const int layer) const { return access()->implementation->mergedModules(layer); } + + /// set the coordinate offset in theta + inline void setOffsetTheta(double offset) const { access()->implementation->setOffsetTheta(offset); } + + /// set the grid size in theta + inline void setGridSizeTheta(double cellSize) const { access()->implementation->setGridSizeTheta(cellSize); } + + /// access the field name used for theta + inline const std::string& fieldNameTheta() const { return access()->implementation->fieldNameTheta(); } + + /// access the field name used for module + inline const std::string& fieldNameModule() const { return access()->implementation->fieldNameModule(); } + + /// access the field name used for layer + inline const std::string& fieldNameLayer() const { return access()->implementation->fieldNameLayer(); } + + + /** \brief Returns a std::vector of the cellDimensions of the given cell ID + in natural order of dimensions (nModules, dTheta). Not implemented yet. + + Returns a std::vector of the cellDimensions of the given cell ID + \param cellID + \return std::vector size 2: + -# size in module + -# size in theta + Not implemented yet. + */ + inline std::vector cellDimensions(const CellID& /*id*/) const { + //return {access()->implementation->gridSizePhi(), access()->implementation->gridSizeTheta()}; + // not implemented + throw std::runtime_error( Form("Function %s not implemented", __PRETTY_FUNCTION__) ); + + return {0., 0.}; + } +}; + +} /* End namespace dd4hep */ +#endif // DD4HEP_DDCORE_GRIDMODULETHETAMERGED_H diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h new file mode 100644 index 00000000..796ba400 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h @@ -0,0 +1,117 @@ +#ifndef DETSEGMENTATION_FCCSWGridModuleThetaMerged_k4geo_H +#define DETSEGMENTATION_FCCSWGridModuleThetaMerged_k4geo_H + +// FCCSW +#include "detectorSegmentations/GridTheta_k4geo.h" +#include "DD4hep/VolumeManager.h" + +/** FCCSWGridModuleThetaMerged_k4geo Detector/DetSegmentation/DetSegmentation/FCCSWGridModuleThetaMerged_k4geo.h FCCSWGridModuleThetaMerged_k4geo.h + * + * Segmentation in theta and module. + * Based on GridTheta, merges modules and theta cells based on layer number + * + */ + +namespace dd4hep { +namespace DDSegmentation { +class FCCSWGridModuleThetaMerged_k4geo : public GridTheta_k4geo { +public: + /// default constructor using an arbitrary type + FCCSWGridModuleThetaMerged_k4geo(const std::string& aCellEncoding); + /// Default constructor used by derived classes passing an existing decoder + FCCSWGridModuleThetaMerged_k4geo(const BitFieldCoder* decoder); + + /// destructor + virtual ~FCCSWGridModuleThetaMerged_k4geo() = default; + + /// read n(modules) from detector metadata + void GetNModulesFromGeom(); + + /// read n(layers) from detector metadata + void GetNLayersFromGeom(); + + /** Determine the local position based on the cell ID. + * @param[in] aCellId ID of a cell. + * return Position (relative to R, phi of Geant4 volume it belongs to, scaled for R=1). + */ + virtual Vector3D position(const CellID& aCellID) const; + /** Determine the cell ID based on the position. + * @param[in] aLocalPosition (not used). + * @param[in] aGlobalPosition + * @param[in] aVolumeId ID of the Geant4 volume + * return Cell ID. + */ + virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, + const VolumeID& aVolumeID) const; + /** Determine the azimuthal angle (relative to the G4 volume) based on the cell ID. + * @param[in] aCellId ID of a cell. + * return Phi. + */ + double phi(const CellID& aCellID) const; + /** Determine the polar angle (relative to the G4 volume) based on the cell ID. + * @param[in] aCellId ID of a cell. + * return Theta. + */ + double theta(const CellID& aCellID) const; + /** Determine the radius based on the cell ID. + * @param[in] aCellId ID of a cell. + * return Radius. + */ + // double radius(const CellID& aCellID) const; + /** Get the number of merged cells in theta for given layer + * @param[in] layer + * return The number of merged cells in theta + */ + inline int mergedThetaCells(const int layer) const { + if (layer < (int) m_mergedCellsTheta.size()) + return m_mergedCellsTheta[layer]; + else + return 1; + } + /** Get the number of merged modules (inclined in phi) + * @param[in] layer + * return The number of merged modules + */ + inline int mergedModules(const int layer) const { + if (layer < (int) m_mergedModules.size()) + return m_mergedModules[layer]; + else + return 1; + } + /** Get the total number of modules of detector + * return The number of modules (as it was set by the user in the xml file..) + */ + inline int nModules() const { return m_nModules; } + /** Get the number of layers + * return The number of layers + */ + inline int nLayers() const { return m_nLayers; } + /** Get the field name used for the layer + * return The field name for the layer. + */ + inline const std::string& fieldNameLayer() const { return m_layerID; } + /** Get the field name used for the module number + * return The field name for the module. + */ + inline const std::string& fieldNameModule() const { return m_moduleID; } + +protected: + /// the field name used for layer + std::string m_layerID; + /// the field name used for the read-out module (can differ from module due to merging) + std::string m_moduleID; + /// vector of number of cells to be merged along theta for each layer (typically between 1 and 4) + std::vector m_mergedCellsTheta; + /// vector of number of modules to be merged for each layer (typically 1 or 2) + std::vector m_mergedModules; + + /// number of modules (or, equivalently, the deltaPhi between adjacent modules) + int m_nModules; + + /// number of layers (from the geometry) + int m_nLayers; + +}; +} +} +#endif /* DETSEGMENTATION_FCCSWGridModuleThetaMerged_k4geo_H */ diff --git a/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h index 038fdf06..84236c9d 100644 --- a/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h @@ -69,7 +69,7 @@ class GridTheta_k4geo : public Segmentation { inline void setFieldNameTheta(const std::string& fieldName) { m_thetaID = fieldName; } /// calculates the Cartesian position from cylindrical coordinates (r, phi, theta) inline Vector3D positionFromRThetaPhi(double ar, double atheta, double aphi) const { - return Vector3D(ar * std::cos(aphi), ar * std::sin(aphi), ar / std::tan(atheta)); + return Vector3D(ar * std::cos(aphi), ar * std::sin(aphi), ar * std::cos(atheta)/std::sin(atheta)); } /// calculates the theta angle from Cartesian coordinates inline double thetaFromXYZ(const Vector3D& aposition) const { diff --git a/detectorSegmentations/src/FCCSWGridModuleThetaMergedHandle_k4geo.cpp b/detectorSegmentations/src/FCCSWGridModuleThetaMergedHandle_k4geo.cpp new file mode 100644 index 00000000..c70004a6 --- /dev/null +++ b/detectorSegmentations/src/FCCSWGridModuleThetaMergedHandle_k4geo.cpp @@ -0,0 +1,4 @@ +#include "detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h" +#include "DD4hep/detail/Handle.inl" + +DD4HEP_INSTANTIATE_HANDLE_UNNAMED(dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo); diff --git a/detectorSegmentations/src/FCCSWGridModuleThetaMerged_k4geo.cpp b/detectorSegmentations/src/FCCSWGridModuleThetaMerged_k4geo.cpp new file mode 100644 index 00000000..ac54d86a --- /dev/null +++ b/detectorSegmentations/src/FCCSWGridModuleThetaMerged_k4geo.cpp @@ -0,0 +1,156 @@ +#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" + +#include +#include "DD4hep/Detector.h" + +namespace dd4hep { +namespace DDSegmentation { + +/// default constructor using an encoding string +FCCSWGridModuleThetaMerged_k4geo::FCCSWGridModuleThetaMerged_k4geo(const std::string& cellEncoding) : GridTheta_k4geo(cellEncoding) { + // define type and description + _type = "FCCSWGridModuleThetaMerged_k4geo"; + _description = "Module-theta segmentation with per-layer merging along theta and/or module"; + + // register all necessary parameters (additional to those registered in GridTheta_k4geo) + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); + registerIdentifier("identifier_module", "Cell ID identifier for readout module", m_moduleID, "module"); + registerParameter("mergedCells_Theta", "Numbers of merged cells in theta per layer", m_mergedCellsTheta, std::vector()); + registerParameter("mergedModules", "Numbers of merged modules per layer", m_mergedModules, std::vector()); + GetNModulesFromGeom(); + GetNLayersFromGeom(); +} + +FCCSWGridModuleThetaMerged_k4geo::FCCSWGridModuleThetaMerged_k4geo(const BitFieldCoder* decoder) : GridTheta_k4geo(decoder) { + // define type and description + _type = "FCCSWGridModuleThetaMerged_k4geo"; + _description = "Module-theta segmentation with per-layer merging along theta and/or module"; + + // register all necessary parameters (additional to those registered in GridTheta_k4geo) + registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer"); + registerIdentifier("identifier_module", "Cell ID identifier for module", m_moduleID, "module"); + registerParameter("mergedCells_Theta", "Numbers of merged cells in theta per layer", m_mergedCellsTheta, std::vector()); + registerParameter("mergedModules", "Numbers of merged modules per layer", m_mergedModules, std::vector()); + GetNModulesFromGeom(); + GetNLayersFromGeom(); +} + +void FCCSWGridModuleThetaMerged_k4geo::GetNModulesFromGeom() { + dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); + try { + m_nModules = dd4hepgeo->constant("ECalBarrelNumPlanes"); + } + catch(...) { + std::cout << "Number of modules not found in detector metadata, exiting..." << std::endl; + exit(1); + } + std::cout << "Number of modules read from detector metadata and used in readout class: " << m_nModules << std::endl; +} + +void FCCSWGridModuleThetaMerged_k4geo::GetNLayersFromGeom() { + dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); + try { + m_nLayers = dd4hepgeo->constant("ECalBarrelNumLayers"); + } + catch(...) { + std::cout << "Number of layers not found in detector metadata, exiting..." << std::endl; + exit(1); + } + std::cout << "Number of layers read from detector metadata and used in readout class: " << m_nLayers << std::endl; +} + +/// determine the local position based on the cell ID +Vector3D FCCSWGridModuleThetaMerged_k4geo::position(const CellID& cID) const { + + // debug + // std::cout << "cellID: " << cID << std::endl; + + // cannot return for R=0 otherwise will lose phi info, return for R=1 + return positionFromRThetaPhi(1.0, theta(cID), phi(cID)); +} + +/// determine the cell ID based on the global position +CellID FCCSWGridModuleThetaMerged_k4geo::cellID(const Vector3D& /* localPosition */, const Vector3D& globalPosition, + const VolumeID& vID) const { + CellID cID = vID; + + // retrieve layer (since merging depends on layer) + int layer = _decoder->get(vID, m_layerID); + + // retrieve theta + double lTheta = thetaFromXYZ(globalPosition); + + // calculate theta bin with original segmentation + int thetaBin = positionToBin(lTheta, m_gridSizeTheta, m_offsetTheta); + + // adjust theta bin if cells are merged along theta in this layer + // assume that m_mergedCellsTheta[layer]>=1 + thetaBin -= (thetaBin % m_mergedCellsTheta[layer]); + + // set theta field of cellID + _decoder->set(cID, m_thetaID, thetaBin); + + // retrieve module number + int module = _decoder->get(vID, m_moduleID); + + // adjust module number if modules are merged in this layer + // assume that m_mergedModules[layer]>=1 + module -= (module % m_mergedModules[layer]); + + // set module field of cellID + _decoder->set(cID, m_moduleID, module); + + return cID; +} + +/// determine the azimuth based on the cell ID +/// the value returned is the relative shift in phi +/// with respect to the first module in the group of +/// merged ones - which will be then added on top of +/// the phi of the volume containing the first cell +/// by the positioning tool +double FCCSWGridModuleThetaMerged_k4geo::phi(const CellID& cID) const { + + // retrieve layer + int layer = _decoder->get(cID, m_layerID); + + // calculate phi offset due to merging + // assume that m_mergedModules[layer]>=1 + double phi = (m_mergedModules[layer]-1) * M_PI / m_nModules ; + + // debug + // std::cout << "layer: " << layer << std::endl; + // std::cout << "merged modules: " << m_mergedModules[layer] << std::endl; + // std::cout << "phi: " << phi << std::endl; + + return phi; +} + +/// determine the polar angle based on the cell ID and the +/// number of merged theta cells +double FCCSWGridModuleThetaMerged_k4geo::theta(const CellID& cID) const { + + // retrieve layer + int layer = _decoder->get(cID, m_layerID); + + // retrieve theta bin from cellID and determine theta position + CellID thetaValue = _decoder->get(cID, m_thetaID); + double _theta = binToPosition(thetaValue, m_gridSizeTheta, m_offsetTheta); + + // adjust return value if cells are merged along theta in this layer + // shift by (N-1)*half theta grid size + // assume that m_mergedCellsTheta[layer]>=1 + _theta += (m_mergedCellsTheta[layer]-1) * m_gridSizeTheta / 2.0 ; + + // debug + // std::cout << "layer: " << layer << std::endl; + // std::cout << "theta bin: " << thetaValue << std::endl; + // std::cout << "gridSizeTheta, offsetTheta: " << m_gridSizeTheta << " , " << m_offsetTheta << std::endl; + // std::cout << "merged cells: " << m_mergedCellsTheta[layer] << std::endl; + // std::cout << "theta: " << _theta << std::endl; + + return _theta; +} + +} +} diff --git a/detectorSegmentations/src/GridEta_k4geo.cpp b/detectorSegmentations/src/GridEta_k4geo.cpp index 064613a0..fada906e 100644 --- a/detectorSegmentations/src/GridEta_k4geo.cpp +++ b/detectorSegmentations/src/GridEta_k4geo.cpp @@ -11,7 +11,7 @@ GridEta_k4geo::GridEta_k4geo(const std::string& cellEncoding) : Segmentation(cel _description = "Eeta segmentation in the global coordinates"; // register all necessary parameters - registerParameter("grid_size_eta", "Cell size in Eta", m_gridSizeEta, 1., SegmentationParameter::LengthUnit); + registerParameter("grid_size_eta", "Cell size in eta", m_gridSizeEta, 1., SegmentationParameter::LengthUnit); registerParameter("offset_eta", "Angular offset in eta", m_offsetEta, 0., SegmentationParameter::AngleUnit, true); registerIdentifier("identifier_eta", "Cell ID identifier for eta", m_etaID, "eta"); } @@ -19,10 +19,10 @@ GridEta_k4geo::GridEta_k4geo(const std::string& cellEncoding) : Segmentation(cel GridEta_k4geo::GridEta_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { // define type and description _type = "GridEta_k4geo"; - _description = "Eeta segmentation in the global coordinates"; + _description = "Eta segmentation in the global coordinates"; // register all necessary parameters - registerParameter("grid_size_eta", "Cell size in Eta", m_gridSizeEta, 1., SegmentationParameter::LengthUnit); + registerParameter("grid_size_eta", "Cell size in eta", m_gridSizeEta, 1., SegmentationParameter::LengthUnit); registerParameter("offset_eta", "Angular offset in eta", m_offsetEta, 0., SegmentationParameter::AngleUnit, true); registerIdentifier("identifier_eta", "Cell ID identifier for eta", m_etaID, "eta"); } @@ -46,7 +46,7 @@ CellID GridEta_k4geo::cellID(const Vector3D& /* localPosition */, const Vector3D // return binToPosition(etaValue, m_gridSizeEta, m_offsetEta); //} -/// determine the polar angle theta based on the cell ID +/// determine the pseudorapidity based on the cell ID double GridEta_k4geo::eta(const CellID& cID) const { CellID etaValue = _decoder->get(cID, m_etaID); return binToPosition(etaValue, m_gridSizeEta, m_offsetEta); diff --git a/detectorSegmentations/src/plugins/SegmentationFactories.cpp b/detectorSegmentations/src/plugins/SegmentationFactories.cpp index 5591f66c..4dd15369 100644 --- a/detectorSegmentations/src/plugins/SegmentationFactories.cpp +++ b/detectorSegmentations/src/plugins/SegmentationFactories.cpp @@ -17,6 +17,9 @@ DECLARE_SEGMENTATION(GridTheta_k4geo, create_segmentation) +#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" +DECLARE_SEGMENTATION(FCCSWGridModuleThetaMerged_k4geo, create_segmentation) + #include "detectorSegmentations/FCCSWGridPhiEta_k4geo.h" DECLARE_SEGMENTATION(FCCSWGridPhiEta_k4geo, create_segmentation) @@ -26,3 +29,4 @@ DECLARE_SEGMENTATION(GridRPhiEta_k4geo, create_segmentation) +