diff --git a/athena.xml b/athena.xml
index 128cf4c7f81467a126ff359c638fc47f7103a7c7..a896d823c484a59a49fa46a2b5ebdd8742273fdc 100644
--- a/athena.xml
+++ b/athena.xml
@@ -190,7 +190,7 @@
<include ref="compact/ce_mrich.xml"/>
<include ref="compact/tof_endcap.xml"/>
<include ref="compact/forward_trd.xml"/>
- <include ref="compact/forward_rich.xml"/>
+ <include ref="compact/gaseous_rich.xml"/>
<include ref="ip6/B0_tracker.xml"/>
<include ref="ip6/far_forward_offM_tracker.xml"/>
diff --git a/compact/forward_rich.xml b/compact/forward_rich.xml
deleted file mode 100644
index f29d6bf682ca37c36a8fee1edf644471ebc064fc..0000000000000000000000000000000000000000
--- a/compact/forward_rich.xml
+++ /dev/null
@@ -1,68 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<lccdd>
- <define>
- <constant name="ForwardRICH_zmin" value="BarrelTracking_length/2.0 + ForwardTracking_length "/>
- <constant name="ForwardRICH_rmin" value="ForwardPID_rmin1"/>
- <constant name="ForwardRICH_rmax0" value="TrackerBarrel_rmax"/>
- <constant name="ForwardRICH_rmax1" value="EcalBarrel_rmin"/>
- <constant name="ForwardRICH_rmax2" value="Solenoid_rmin-6*cm"/>
- <constant name="ForwardRICHDepth" value="0.9*m"/>
- </define>
-
- <detectors>
- <detector
- id="ForwardRICH_ID"
- name="ForwardRICH"
- type="refdet_ForwardRICH"
- readout="ForwardRICHHits"
- vis="RICHVis">
- <dimensions
- z0="ForwardRICH_zmin"
- snout_length="ForwardRICH_length - ForwardRICHDepth"
- length="ForwardRICH_length"
- rmin="ForwardRICH_rmin"
- rmax0="ForwardRICH_rmax0"
- rmax1="ForwardRICH_rmax1"
- rmax2="ForwardRICH_rmax2"/>
- <radiator material="N2cherenkov" />
- <tank
- length="ForwardRICHDepth"
- gas="N2cherenkov"
- vis="GreenVis"
- rmin="ForwardRICH_rmin"
- rmax1="ForwardRICH_rmax1"
- rmax2="ForwardRICH_rmax2" />
- <comment> What are the following MCP-PMT parameters?</comment>
- <mcppmt
- z0="0.0*cm"
- rmin="ForwardRICH_rmax1"
- rmax="ForwardRICH_rmax2"
- rtol="1.0*cm"
- vis="BlueVis"
- module_size="10*cm"
- module_gap="0.2*cm"
- thickness="1.0*cm"
- material="Quartz" />
- <mirror
- z0="ForwardRICH_length - 40*cm"
- thickness="1*mm"
- material="PyrexGlass"
- vis="GrayVis">
- <slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="0*degree" />
- <slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="60*degree" />
- <slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="120*degree" />
- <slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="180*degree" />
- <slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="240*degree" />
- <slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="300*degree" />
- </mirror>
- </detector>
- </detectors>
-
- <readouts>
- <readout name="ForwardRICHHits">
- <segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" />
- <id>system:8,layer:4,piece:4,module:14,x:32:-16,y:-16</id>
- </readout>
- </readouts>
-
-</lccdd>
diff --git a/compact/gaseous_rich.xml b/compact/gaseous_rich.xml
new file mode 100644
index 0000000000000000000000000000000000000000..30e76efe985ceb27112097df28173a9b6e117da6
--- /dev/null
+++ b/compact/gaseous_rich.xml
@@ -0,0 +1,73 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+ <define>
+ <constant name="ForwardRICH_zmin" value="BarrelTracking_length/2.0 + ForwardTracking_length "/>
+ <constant name="ForwardRICH_rmin" value="ForwardPID_rmin1"/>
+ <constant name="ForwardRICH_rmax0" value="TrackerBarrel_rmax"/>
+ <constant name="ForwardRICH_rmax1" value="EcalBarrel_rmin"/>
+ <constant name="ForwardRICH_rmax2" value="Solenoid_rmin-6*cm"/>
+ <constant name="ForwardRICHDepth" value="0.9*m"/>
+ <constant name="ForwardRICH_sensor_rmin" value="ForwardRICH_rmax1 - 40.*cm"/>
+ <constant name="ForwardRICH_sensor_rmax" value="ForwardRICH_rmax2"/>
+ </define>
+
+ <detectors>
+ <detector
+ id="ForwardRICH_ID"
+ name="GaseousRICH"
+ type="athena_GaseousRICH"
+ readout="ForwardRICHHits"
+ gas="N2cherenkov"
+ vis="RICHVis">
+ <dimensions
+ z0="ForwardRICH_zmin"
+ snout_length="ForwardRICH_length - ForwardRICHDepth"
+ length="ForwardRICH_length"
+ rmin="ForwardRICH_rmin"
+ rmax0="ForwardRICH_rmax0"
+ rmax1="ForwardRICH_rmax1"
+ rmax2="ForwardRICH_rmax2"/>
+ <mirrors thickness="1*mm" material="PyrexGlass" vis="GrayVis">
+ <position z="ForwardRICH_length - 13*cm"/>
+ <mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
+ rotz="0*degree" roty="-15*degree" rotx="0"/>
+ <mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
+ rotz="60*degree" roty="-15*degree" rotx="0"/>
+ <mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
+ rotz="120*degree" roty="-15*degree" rotx="0"/>
+ <mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
+ rotz="180*degree" roty="-15*degree" rotx="0"/>
+ <mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
+ rotz="240*degree" roty="-15*degree" rotx="0"/>
+ <mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
+ rotz="300*degree" roty="-15*degree" rotx="0"/>
+ </mirrors>
+ <sensors>
+ <module sx="10*cm" sy="10*cm" sz="1.0*cm" gap="0.2*cm" material="Quartz" vis="AnlGold"/>
+ <comment>A thin optical material to accept optical photon in simulation</comment>
+ <optical material="AirOptical" thickness="0.1*mm"/>
+ <position z="8.*cm"/>
+ <sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
+ rotz="0*degree" roty="-15*degree" rotx="0"/>
+ <sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
+ rotz="60*degree" roty="-15*degree" rotx="0"/>
+ <sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
+ rotz="120*degree" roty="-15*degree" rotx="0"/>
+ <sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
+ rotz="180*degree" roty="-15*degree" rotx="0"/>
+ <sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
+ rotz="240*degree" roty="-15*degree" rotx="0"/>
+ <sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
+ rotz="300*degree" roty="-15*degree" rotx="0"/>
+ </sensors>
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="ForwardRICHHits">
+ <segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" />
+ <id>system:8,sector:8,module:16,x:32:-16,y:-16</id>
+ </readout>
+ </readouts>
+
+</lccdd>
diff --git a/src/ForwardRICH_geo.cpp b/src/ForwardRICH_geo.cpp
deleted file mode 100644
index 0c8f73b06c8419307c5d17d3670871cdb152b7af..0000000000000000000000000000000000000000
--- a/src/ForwardRICH_geo.cpp
+++ /dev/null
@@ -1,192 +0,0 @@
-//==========================================================================
-// Forward Ring Imaging Cherenkov Detector
-//--------------------------------------------------------------------------
-//
-// Author: C. Peng (ANL)
-// Date: 09/30/2020
-//
-//==========================================================================
-
-#include <XML/Helper.h>
-#include "TMath.h"
-#include "TString.h"
-#include "GeometryHelpers.h"
-#include "Math/Point2D.h"
-#include "DDRec/Surface.h"
-#include "DDRec/DetectorData.h"
-#include "DD4hep/OpticalSurfaces.h"
-#include "DD4hep/DetFactoryHelper.h"
-#include "DD4hep/Printout.h"
-
-using namespace std;
-using namespace dd4hep;
-using namespace dd4hep::rec;
-
-
-// create the detector
-static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
-{
- xml::DetElement detElem = handle;
-
- std::string detName = detElem.nameStr();
- int detID = detElem.id();
-
- DetElement det(detName, detID);
- xml::Component dims = detElem.dimensions();
- xml::Component rads = detElem.child(_Unicode(radiator));
- xml::Component mir = detElem.child(_Unicode(mirror));
- xml::Component mcp = detElem.child(_Unicode(mcppmt));
- xml::Component tank = detElem.child(_Unicode(tank));
-
- // dimensions
- double z0 = dims.z0();
- double length = dims.length();
- double rmin = dims.rmin();
- double rmax0 = dims.attr<double>(_Unicode(rmax0));
- double rmax1 = dims.attr<double>(_Unicode(rmax1));
- double rmax2 = dims.attr<double>(_Unicode(rmax2));
- double snout_length = dims.attr<double>(_Unicode(snout_length));
-
- // mirror setting
- auto mThick = mir.thickness();
- auto mirZ = mir.attr<double>(_Unicode(z0));
-
- // mcppmt setting
- auto pRmin = mcp.rmin();
- auto pRmax = mcp.rmax();
- auto pThick = mcp.thickness();
- auto pSize = mcp.attr<double>(_Unicode(module_size));
- auto pGap = mcp.attr<double>(_Unicode(module_gap));
- auto pTol = mcp.attr<double>(_Unicode(rtol));
- auto pZ = mcp.attr<double>(_Unicode(z0));
-
- // tank parameters
- double tank_length = length - snout_length;
-
- // materials
- auto mirMat = desc.material(mir.materialStr());
- auto gasMat = desc.material(rads.materialStr());
- auto mcpMat = desc.material(mcp.materialStr());
-
- double front_offset = snout_length+tank_length/2.0;
-
- // constants
- auto richCenterAngle = std::atan((rmin + (rmax1 - rmin)/2.)/(front_offset+mirZ));
- //std::cout << richCenterAngle*180./M_PI << std::endl;
-
- // an envelope for the detector
- // use a complicated shape to avoid conflict with the other parts
- // cone for radiator and the first set of mirrors
- double halfLength = length/2.;
- Cone env1(snout_length/2.0, rmin, rmax0, rmin, rmax1);
- // envelope for detection plane
- // Cone env2(halfLength - pZ/2., rmin, pRmax, rmin, rmax2);
- Tube env2(rmin, pRmax + pTol + pGap + 1.0*cm, tank_length/2., 0., 2*M_PI);
-
- UnionSolid envShape(env2, env1, Position(0., 0., -tank_length/2.-snout_length/2));
-
- Volume envVol(detName + "_envelope", envShape, gasMat);
- envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
-
- // ---------------
- // spherical mirrors inside it
- int ilayer = 1;
-
- // optical surface
- OpticalSurfaceManager surfMgr = desc.surfaceManager();
- OpticalSurface mirSurf = surfMgr.opticalSurface("MirrorOpticalSurface");
- // mirror slices
- int imod = 1;
- for (xml::Collection_t sl(mir, _Unicode(slice)); sl; ++sl, ++imod) {
- auto focus = sl.attr<double>(_Unicode(focus));
- auto wphi = sl.attr<double>(_Unicode(phiw));
- auto rotZ = sl.attr<double>(_Unicode(rotz));
- auto mRmin = sl.attr<double>(_Unicode(rmin));
- auto mRmax = sl.attr<double>(_Unicode(rmax));
- double curve = 0.;
- if (sl.hasAttr(_Unicode(curve))) {
- curve = sl.attr<double>(_Unicode(curve));
- }
- // geometry of mirror slice
- PlacedVolume mirPV;
- Volume mirVol(Form("mirror_v_dummy%d", imod));
- mirVol.setMaterial(mirMat);
- mirVol.setVisAttributes(desc.visAttributes(mir.visStr()));
- // spherical mirror
- if (curve > 0.) {
- // somehow geant4 does not support -wphi/2. to wphi/2., so additonal rotation in Z
- double mTheta1 = std::asin(mRmin/curve);
- double mTheta2 = std::asin(mRmax/curve);
- double rotY = -std::asin(focus/curve);
- mirVol.setSolid(Sphere(curve, curve + mThick, mTheta1*1.01, mTheta2*0.99, 0., wphi));
- // action is in a reverse order
- Transform3D tr = Translation3D(0., 0., mirZ - front_offset) // move for z position
- * RotationZ(rotZ) // rotate phi angle
- * RotationY(rotY) // rotate for focus point
- * RotationX(180*degree)
- * Translation3D(0., 0., -curve) // move spherical shell to origin
- * RotationZ(-wphi/2.); // center phi angle to 0. (-wphi/2., wphi/2.)
- mirPV = envVol.placeVolume(mirVol, tr);
- // plane mirror
- } else {
- mirVol.setSolid(Tube(mRmin, mRmax, mThick/2.0, 0., wphi));
- Transform3D tr = Translation3D(0., 0., mirZ - front_offset) // move for z position
- * RotationZ(rotZ) // rotate phi angle
- * RotationZ(-wphi/2.); // center phi angle to 0. (-wphi/2., wphi/2.)
- mirPV = envVol.placeVolume(mirVol, tr);
- }
- mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", imod);
- DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
- mirDE.setPlacement(mirPV);
- SkinSurface mirSurfBorder(desc, mirDE, Form("RICHmirror%d", imod), mirSurf, mirVol);
- mirSurfBorder.isValid();
- }
- ilayer++;
-
- // ---------------
- // photo-detector unit
- // Fill the photo-detection plane with square shape MCP-PMTs
- Box mcpShape1(pSize/2.0, pSize/2.0, pThick/2.0);
- Volume mcpVol1("mcppmt_v_material", mcpShape1, mcpMat);
-
- // a thin layer of cherenkov gas for accepting optical photons
- Box mcpShape(pSize/2.0, pSize/2.0, pThick/2.0 + 0.1*mm);
- Volume mcpVol("mcppmt_v", mcpShape, gasMat);
- mcpVol.placeVolume(mcpVol1, Position(0., 0., -0.1*mm));
-
- mcpVol.setVisAttributes(desc.visAttributes(mcp.visStr()));
- sens.setType("photoncounter");
- mcpVol.setSensitiveDetector(sens);
-
- // photo-detector plane envelope
- for (size_t ipd = 0; ipd < 6; ++ipd) {
- double phmin = -M_PI/6.5; // added 0.5 to make it smaller
- double phmax = M_PI/6.5;
- Tube pdEnvShape(pRmin - pTol - pGap, pRmax + pTol + pGap, pThick/2.0 + 0.1*cm, phmin, phmax);
- Volume pdVol("pd_envelope", pdEnvShape, desc.material("AirOptical"));
- auto points = ref::utils::fillSquares({0., 0.}, pSize + pGap, pRmin + pTol + pGap, pRmax + pTol + pGap, phmin, phmax);
- for (size_t i = 0; i < points.size(); ++i) {
- auto pt = points[i];
- auto mcpPV = pdVol.placeVolume(mcpVol, Position(pt.x(), pt.y(), 0.));
- mcpPV.addPhysVolID("layer", ilayer).addPhysVolID("module", i + 1);
- DetElement mcpDE(det, Form("MCPPMT_DE%d_%d", ipd + 1, i + 1), i + 1);
- mcpDE.setPlacement(mcpPV);
- }
- Transform3D tr = Translation3D(0., 0., -front_offset + pZ + pThick/2.0) // move for z position
- * RotationZ(ipd*M_PI/3.) // rotate phi angle
- * RotationY(-richCenterAngle); // rotate to perpendicular position
- auto pdPV = envVol.placeVolume(pdVol, tr);
- pdPV.addPhysVolID("layer", ilayer).addPhysVolID("piece", ipd + 1);
- }
- Volume motherVol = desc.pickMotherVolume(det);
- PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0 + front_offset));
- envPV.addPhysVolID("system", detID);
- det.setPlacement(envPV);
-
- return det;
-}
-//@}
-
-// clang-format off
-DECLARE_DETELEMENT(refdet_ForwardRICH, createDetector)
-
diff --git a/src/GaseousRICH_geo.cpp b/src/GaseousRICH_geo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..38a66b8de76ede7147e64147b088712fcdf96a0e
--- /dev/null
+++ b/src/GaseousRICH_geo.cpp
@@ -0,0 +1,215 @@
+//==========================================================================
+// Gaseous Ring Imaging Cherenkov Detector
+//--------------------------------------------------------------------------
+//
+// Author: C. Peng (ANL)
+// Date: 09/30/2020
+//
+//==========================================================================
+
+#include <XML/Helper.h>
+#include "TMath.h"
+#include "TString.h"
+#include "GeometryHelpers.h"
+#include "Math/Point2D.h"
+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+// headers
+void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset);
+void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset);
+void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens);
+
+// helper function to get x, y, z if defined in a xml component
+template<class XmlComp>
+Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
+{
+ Position pos(0., 0., 0.);
+ if (comp.hasChild(name)) {
+ auto child = comp.child(name);
+ pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
+ pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
+ pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
+ }
+ return pos;
+}
+
+// create the detector
+static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
+{
+ xml::DetElement detElem = handle;
+ std::string detName = detElem.nameStr();
+ int detID = detElem.id();
+
+ DetElement det(detName, detID);
+ xml::Component dims = detElem.dimensions();
+
+ // build a big envelope for all the components, filled with optical material to ensure transport of optical photons
+ double z0 = dims.z0();
+ double length = dims.length();
+ double rmin = dims.rmin();
+ double rmax0 = dims.attr<double>(_Unicode(rmax0));
+ double rmax1 = dims.attr<double>(_Unicode(rmax1));
+ double rmax2 = dims.attr<double>(_Unicode(rmax2));
+ double snout_length = dims.attr<double>(_Unicode(snout_length));
+ // fill envelope with radiator materials (need optical property for optical photons)
+ auto gasMat = desc.material(dd4hep::getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
+
+ Cone snout(snout_length/2.0, rmin, rmax0, rmin, rmax1);
+ Tube tank(rmin, rmax2, (length - snout_length)/2., 0., 2*M_PI);
+ // shift the snout to the left side of tank
+ UnionSolid envShape(tank, snout, Position(0., 0., -length/2.));
+ // some reference points
+ auto snout_front = Position(0., 0., -(length + snout_length)/2.);
+ // auto snout_end = Position(0., 0., -(length - snout_length)/2.); // tank_front
+ // auto tank_end = Position(0., 0., (length - snout_length)/2.);
+
+ Volume envVol(detName + "_envelope", envShape, gasMat);
+ envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
+
+ // sensitive detector type
+ sens.setType("photoncounter");
+
+ // @TODO: place a radiator
+ // build_radiator(desc, envVol, detElement.child(_Unicode(radiator)), snout_front);
+
+ // place mirrors
+ build_mirrors(desc, det, envVol, detElem.child(_Unicode(mirrors)), snout_front);
+
+ // place photo-sensors
+ build_sensors(desc, envVol, detElem.child(_Unicode(sensors)), snout_front, sens);
+
+ Volume motherVol = desc.pickMotherVolume(det);
+ PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0) - snout_front);
+ envPV.addPhysVolID("system", detID);
+ det.setPlacement(envPV);
+
+ return det;
+}
+
+// @TODO: implement a radiator, now the envelope serves as the radiator
+void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset)
+{
+ // place holder
+}
+
+// place mirrors
+void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset)
+{
+ double thickness = dd4hep::getAttrOrDefault<double>(plm, _Unicode(dz), 1.*dd4hep::mm);
+ auto mat = desc.material(plm.attr<std::string>(_Unicode(material)));
+ auto vis = desc.visAttributes(plm.attr<std::string>(_Unicode(vis)));
+
+ // optical surface
+ OpticalSurfaceManager surfMgr = desc.surfaceManager();
+ auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(plm, _Unicode(surface), "MirrorOpticalSurface"));
+
+ // placements
+ auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
+ auto grot = get_xml_xyz(plm, _Unicode(position));
+ int imod = 1;
+ for (xml::Collection_t mir(plm, _Unicode(mirror)); mir; ++mir, ++imod) {
+ double rmin = mir.attr<double>(_Unicode(rmin));
+ double rmax = mir.attr<double>(_Unicode(rmax));
+ double phiw = dd4hep::getAttrOrDefault<double>(mir, _Unicode(phiw), 2.*M_PI);
+ double rotz = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotz), 0.);
+ double roty = dd4hep::getAttrOrDefault<double>(mir, _Unicode(roty), 0.);
+ double rotx = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotx), 0.);
+
+ Volume vol(Form("mirror_v_dummy%d", imod));
+ vol.setMaterial(mat);
+ vol.setVisAttributes(vis);
+ // mirror curvature
+ double curve = dd4hep::getAttrOrDefault<double>(mir, _Unicode(curve), 0.);
+ // spherical mirror
+ if (curve > 0.) {
+ double th1 = std::asin(rmin/curve);
+ double th2 = std::asin(rmax/curve);
+ vol.setSolid(Sphere(curve, curve + thickness, th1, th2, 0., phiw));
+ // plane mirror
+ } else {
+ vol.setSolid(Tube(rmin, rmax, thickness/2., 0., phiw));
+ }
+ // transforms are in a reverse order
+ Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
+ * RotationZYX(grot.z(), grot.y(), grot.x())
+ * RotationZ(rotz) // rotation of the piece
+ * RotationY(roty) // rotation of the piece
+ * RotationX(rotx) // rotation of the piece
+ * Translation3D(0., 0., -curve) // move spherical shell to origin (no move for planes)
+ * RotationZ(-phiw/2.); // center phi angle to 0. (-phiw/2., phiw/2.)
+ auto pv = env.placeVolume(vol, tr);
+ DetElement de(sdet, Form("mirror_de%d", imod), imod);
+ de.setPlacement(pv);
+ SkinSurface skin(desc, de, Form("mirror_optical_surface%d", imod), surf, vol);
+ skin.isValid();
+ }
+}
+
+// place photo-sensors
+void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens)
+{
+ // build sensor unit geometry
+ auto mod = plm.child(_Unicode(module));
+ double sx = mod.attr<double>(_Unicode(sx));
+ double sy = mod.attr<double>(_Unicode(sy));
+ double sz = mod.attr<double>(_Unicode(sz));
+ double gap = mod.attr<double>(_Unicode(gap));
+ auto mat = desc.material(mod.attr<std::string>(_Unicode(material)));
+ auto vis = desc.visAttributes(mod.attr<std::string>(_Unicode(vis)));
+
+ Box sensor(sx/2., sy/2., sz/2.);
+ Volume svol("sensor_v", sensor, mat);
+ svol.setVisAttributes(vis);
+
+ // a thin layer of cherenkov gas for accepting optical photons
+ auto opt = plm.child(_Unicode(optical));
+ double opthick = opt.attr<double>(_Unicode(thickness));
+ auto opmat = desc.material(opt.attr<std::string>(_Unicode(material)));
+
+ Box opshape(sx/2., sy/2., sz/2. + opthick/2.);
+ Volume opvol("sensor_v_optical", opshape, opmat);
+ opvol.placeVolume(svol, Position(0., 0., 0.));
+ opvol.setSensitiveDetector(sens);
+
+ // photo-detector plane envelope
+ auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
+ auto grot = get_xml_xyz(plm, _Unicode(position));
+ int isec = 1;
+ for (xml::Collection_t sec(plm, _Unicode(sector)); sec; ++sec, ++isec) {
+ double rmin = sec.attr<double>(_Unicode(rmin));
+ double rmax = sec.attr<double>(_Unicode(rmax));
+ double phiw = dd4hep::getAttrOrDefault<double>(sec, _Unicode(phiw), 2.*M_PI);
+ double rotz = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotz), 0.);
+ double roty = dd4hep::getAttrOrDefault<double>(sec, _Unicode(roty), 0.);
+ double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
+
+ // fill sensors to the piece
+ auto points = ref::utils::fillRectangles({0., 0.}, sx + gap, sy + gap, rmin - gap, rmax + gap, -phiw/2., phiw/2.);
+ int imod = 1;
+ for (auto &p : points) {
+ // transofrms are in a reversed order
+ Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
+ * RotationZYX(grot.z(), grot.y(), grot.x())
+ * RotationZ(rotz) // rotation of the sector
+ * RotationY(roty) // rotation of the sector
+ * RotationX(rotx) // rotation of the sector
+ * Translation3D(p.x(), p.y(), 0.); // place modules in each sector
+ auto pv = env.placeVolume(opvol, tr);
+ pv.addPhysVolID("sector", isec).addPhysVolID("module", imod++);
+ }
+ }
+}
+
+//@}
+
+// clang-format off
+DECLARE_DETELEMENT(athena_GaseousRICH, createDetector)
+
diff --git a/src/GeometryHelpers.cpp b/src/GeometryHelpers.cpp
index 201b88d6d41da7d653eef8825b3af4f0cb92e83d..c37d30bd175bfa1e1334b841fb7eb2c0115b4084 100644
--- a/src/GeometryHelpers.cpp
+++ b/src/GeometryHelpers.cpp
@@ -6,7 +6,7 @@ namespace ref::utils {
typedef ROOT::Math::XYPoint Point;
// check if a square in a ring
-inline bool in_ring(const Point &pt, double side, double rmin, double rmax, double phmin, double phmax)
+inline bool in_ring(const Point &pt, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{
if (pt.r() > rmax || pt.r() < rmin) {
return false;
@@ -14,10 +14,10 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax, doub
// check four corners
std::vector<Point> pts {
- Point(pt.x() - side/2., pt.y() - side/2.),
- Point(pt.x() - side/2., pt.y() + side/2.),
- Point(pt.x() + side/2., pt.y() - side/2.),
- Point(pt.x() + side/2., pt.y() + side/2.),
+ Point(pt.x() - sx/2., pt.y() - sy/2.),
+ Point(pt.x() - sx/2., pt.y() + sy/2.),
+ Point(pt.x() + sx/2., pt.y() - sy/2.),
+ Point(pt.x() + sx/2., pt.y() + sy/2.),
};
for (auto &p : pts) {
if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
@@ -28,11 +28,11 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax, doub
}
// check if a square is overlapped with the others
-inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts)
+inline bool overlap(const Point &pt, double sx, double sy, const std::vector<Point> &pts)
{
for (auto &p : pts) {
- auto pn = (p - pt)/side;
- if ((std::abs(pn.x()) < 1. - 1e-6) && (std::abs(pn.y()) < 1. - 1e-6)) {
+ auto pn = p - pt;
+ if ((std::abs(pn.x()) < (1. - 1e-6)*sx) && (std::abs(pn.y()) < (1. - 1e-6)*sy)) {
return true;
}
}
@@ -40,25 +40,24 @@ inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts)
}
// a helper function to recursively fill square in a ring
-void add_square(Point p, std::vector<Point> &res, double lside, double rmin, double rmax,
- double phmin, double phmax)
+void add_rectangle(Point p, std::vector<Point> &res, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{
// outside of the ring or overlapping
- if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, res)) {
+ if (!in_ring(p, sx, sy, rmin, rmax, phmin, phmax) || overlap(p, sx, sy, res)) {
return;
}
res.emplace_back(p);
// check adjacent squares
- add_square(Point(p.x() + lside, p.y()), res, lside, rmin, rmax, phmin, phmax);
- add_square(Point(p.x() - lside, p.y()), res, lside, rmin, rmax, phmin, phmax);
- add_square(Point(p.x(), p.y() + lside), res, lside, rmin, rmax, phmin, phmax);
- add_square(Point(p.x(), p.y() - lside), res, lside, rmin, rmax, phmin, phmax);
+ add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax);
+ add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax);
+ add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax);
+ add_rectangle(Point(p.x(), p.y() - sy), res, sx, sy, rmin, rmax, phmin, phmax);
}
// fill squares
-std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax, double phmin, double phmax)
+std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{
// convert (0, 2pi) to (-pi, pi)
if (phmax > M_PI) {
@@ -67,13 +66,13 @@ std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax
}
// start with a seed square and find one in the ring
// move to center
- ref = ref - Point(int(ref.x()/lside)*lside, int(ref.y()/lside)*lside);
+ ref = ref - Point(int(ref.x()/sx)*sx, int(ref.y()/sy)*sy);
- auto find_seed = [] (const Point &ref, int n, double side, double rmin, double rmax, double phmin, double phmax) {
+ auto find_seed = [] (const Point &ref, int n, double sx, double sy, double rmin, double rmax, double phmin, double phmax) {
for (int ix = -n; ix < n; ++ix) {
for (int iy = -n; iy < n; ++iy) {
- Point pt(ref.x() + ix*side, ref.y() + iy*side);
- if (in_ring(pt, side, rmin, rmax, phmin, phmax)) {
+ Point pt(ref.x() + ix*sx, ref.y() + iy*sy);
+ if (in_ring(pt, sx, sy, rmin, rmax, phmin, phmax)) {
return pt;
}
}
@@ -82,8 +81,8 @@ std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax
};
std::vector<Point> res;
- ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax, phmin, phmax);
- add_square(ref, res, lside, rmin, rmax, phmin, phmax);
+ ref = find_seed(ref, int(rmax/sx) + 2, sx, sy, rmin, rmax, phmin, phmax);
+ add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax);
return res;
}
diff --git a/src/GeometryHelpers.h b/src/GeometryHelpers.h
index 12cf9a814a64ed978f7c8371a1e3a4c8d36c7543..028c474c833809af6fb3c605b7d956fcfe0d9153 100644
--- a/src/GeometryHelpers.h
+++ b/src/GeometryHelpers.h
@@ -7,8 +7,14 @@ namespace ref::utils {
typedef ROOT::Math::XYPoint Point;
+// fill rectangles in a ring
+std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax,
+ double phmin = -M_PI, double phmax = M_PI);
// fill squares in a ring
-std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax,
- double phmin = -M_PI, double phmax = M_PI);
+inline std::vector<Point> fillSquares(Point ref, double size, double rmin, double rmax,
+ double phmin = -M_PI, double phmax = M_PI)
+{
+ return fillRectangles(ref, size, size, rmin, rmax, phmin, phmax);
+}
} // ref::utils
diff --git a/src/HomogeneousCalorimeter_geo.cpp b/src/HomogeneousCalorimeter_geo.cpp
index 65fb78dddbaab2fd75d41094699b481d78d0d75a..a365801dc469d800e55682a2a63db9d24e6ed5e3 100644
--- a/src/HomogeneousCalorimeter_geo.cpp
+++ b/src/HomogeneousCalorimeter_geo.cpp
@@ -275,7 +275,7 @@ static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, Sens
double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
- auto points = ref::utils::fillSquares({0., 0.}, modSize.x(), rmin, rmax, phimin, phimax);
+ auto points = ref::utils::fillRectangles({0., 0.}, modSize.x(), modSize.y(), rmin, rmax, phimin, phimax);
// placement to mother
auto pos = get_xml_xyz(plm, _Unicode(position));
auto rot = get_xml_xyz(plm, _Unicode(rotation));