diff --git a/compact/ce_mrich.xml b/compact/ce_mrich.xml
new file mode 100644
index 0000000000000000000000000000000000000000..ebd1c473ceddbefb04b9b9ea57f544fb7b8fcec2
--- /dev/null
+++ b/compact/ce_mrich.xml
@@ -0,0 +1,20 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+ <define>
+ </define>
+
+ <detectors>
+ <detector id="ce_MRICH_ID" name="ce_MRICH" type="refdet_ce_MRICH" readout="MRICHHits" vis="BlueVis" material="Air">
+ <dimensions rmin="ce_MRICHRMin" rmax="ce_MRICHRMax" length="ce_MRICHLength" z="ce_MRICHZMin"/>
+ <modules material="Quartz" thickness="10*cm" width="10*cm" gap="1*cm" />
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="MRICHHits">
+ <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/definitions.xml b/compact/definitions.xml
index c7321b9729a1c651422d86b0627449cf6cac8872..a4e31e2fa449bbc9e25eb281305cefcb33623f0c 100644
--- a/compact/definitions.xml
+++ b/compact/definitions.xml
@@ -176,7 +176,7 @@
Forwardtracking ID: 120
Forward RICH ID: 121
- Unused IDs: 122-129
+ Unused IDs: 123-129
</comment>
<constant name="ForwardTracking_ID" value="120"/>
<constant name="ForwardRICH_ID" value="121"/>
@@ -186,9 +186,11 @@
(130-139) Backward reserved
=====================================
- TBD
+ Modular RICH ID: 130
+ Unused IDs: 131-139
</comment>
+ <constant name="ce_MRICH_ID" value="130"/>
<comment>
=====================================
@@ -455,6 +457,17 @@
</comment>
+ <comment>
+ --------------------------
+ ce_MRICH Parameters
+ --------------------------
+ </comment>
+ <constant name="ce_MRICHRMin" value="15*cm"/>
+ <constant name="ce_MRICHRMax" value="100*cm"/>
+ <constant name="ce_MRICHLength" value="15*cm"/>
+ <constant name="ce_MRICHZMin" value="-SolenoidYokeEndcap_zmin"/>
+
+
</define>
diff --git a/src/ForwardRICH_geo.cpp b/src/ForwardRICH_geo.cpp
index 6520c30aac30492a8950c1590661bce57e089e47..3dc6417ba4562a43e22227d8f4a7e9ec969c3bd1 100644
--- a/src/ForwardRICH_geo.cpp
+++ b/src/ForwardRICH_geo.cpp
@@ -10,6 +10,7 @@
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
+#include "ref_utils.h"
#include "Math/Point2D.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
@@ -21,86 +22,6 @@ using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
-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)
-{
- if (pt.r() > rmax || pt.r() < rmin) {
- return false;
- }
-
- // 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.),
- };
- for (auto &p : pts) {
- if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
- return false;
- }
- }
- return true;
-}
-
-// check if a square is overlapped with the others
-inline bool overlap(const Point &pt, double side, 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)) {
- return true;
- }
- }
- return false;
-}
-
-// 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)
-{
- // outside of the ring or overlapping
- if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, 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);
-}
-
-// fill squares
-std::vector<Point> fill_squares(Point ref, double lside, double rmin, double rmax,
- double phmin = 0., double phmax = 2.*M_PI)
-{
- // 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);
-
- auto find_seed = [] (const Point &ref, int n, double side, 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)) {
- return pt;
- }
- }
- }
- return ref;
- };
-
- 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);
- return res;
-}
// create the detector
static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
@@ -235,7 +156,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
double phmax = M_PI/6.;
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 = fill_squares(Point(0., 0.), pSize + pGap, pRmin - pTol - pGap, pRmax + pTol + pGap, phmin, phmax);
+ 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.));
diff --git a/src/ce_MRICH.cpp b/src/ce_MRICH.cpp
index a4648190c30b5746cff07aba69d8e694cb3da600..64b9a682be81ac03ea64b12f99bab347e69ba366 100644
--- a/src/ce_MRICH.cpp
+++ b/src/ce_MRICH.cpp
@@ -2,22 +2,18 @@
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.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"
+#include "ref_utils.h"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
-typedef ROOT::Math::XYPoint Point;
-
-std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax,
- double phmin = 0., double phmax = 2.*M_PI);
void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, SensitiveDetector &sens);
@@ -30,18 +26,17 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
DetElement det(detName, detID);
xml::Component dims = detElem.dimensions();
- xml::Component rads = detElem.child(_Unicode(radiator));
+ // xml::Component rads = detElem.child(_Unicode(radiator));
- auto RIn = dims.attr<double>(_Unicode(r_in));
- auto ROut = dims.attr<double>(_Unicode(r_out));
- auto SizeZ = dims.attr<double>(_Unicode(size_z));
- auto PosZ = dims.z();
- auto InnerR = dims.attr<double>(_Unicode(inner_r));
+ auto rmin = dims.rmin();
+ auto rmax = dims.rmax();
+ auto length = dims.length();
+ auto z0 = dims.z();
- auto gasMat = desc.material(detElem.materialStr());
+ auto gasMat = desc.material("Air");
// detector envelope
- Tube envShape(RIn, ROut, SizeZ / 2., 0., 2*M_PI);
+ Tube envShape(rmin, rmax, length/2., 0., 2*M_PI);
Volume envVol("ce_MRICH_GVol", envShape, gasMat);
envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
@@ -50,7 +45,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
// place envelope
Volume motherVol = desc.pickMotherVolume(det);
- PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, PosZ));
+ PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0));
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
@@ -62,14 +57,15 @@ void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, Sensit
xml::Component dims = detElem.dimensions();
xml::Component mods = detElem.child(_Unicode(modules));
- auto ROut = dims.attr<double>(_Unicode(r_out));
+ auto rmin = dims.rmin();
+ auto rmax = dims.rmax();
+
auto mThick = mods.attr<double>(_Unicode(thickness));
auto mWidth = mods.attr<double>(_Unicode(width));
auto mGap = mods.attr<double>(_Unicode(gap));
- auto InnerR = dims.attr<double>(_Unicode(inner_r));
auto modMat = desc.material(mods.materialStr());
- auto gasMat = desc.material(detElem.materialStr());
+ auto gasMat = desc.material("Air");
// single module
Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm);
@@ -85,11 +81,8 @@ void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, Sensit
sens.setType("photoncounter");
modVol.setSensitiveDetector(sens);
- double MRICH_R, x_mrich = 0., y_mrich = 0;
- int kmrich = -1;
-
// place modules in the sectors (disk)
- auto points = fillSquares(Point(0., 0.), mWidth + mGap, InnerR, ROut + mGap);
+ auto points = ref::utils::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap);
// determine module direction, always facing z = 0
double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.;
@@ -103,83 +96,6 @@ void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, Sensit
}
}
-// check if a square in a ring
-inline bool in_ring(const Point &pt, double side, double rmin, double rmax, double phmin, double phmax)
-{
- if (pt.r() > rmax || pt.r() < rmin) {
- return false;
- }
-
- // 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.),
- };
- for (auto &p : pts) {
- if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
- return false;
- }
- }
- return true;
-}
-
-// check if a square is overlapped with the others
-inline bool overlap(const Point &pt, double side, 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)) {
- return true;
- }
- }
- return false;
-}
-
-// 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)
-{
- // outside of the ring or overlapping
- if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, 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);
-}
-
-// fill squares
-std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax, double phmin, double phmax)
-{
- // 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);
-
- auto find_seed = [] (const Point &ref, int n, double side, 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)) {
- return pt;
- }
- }
- }
- return ref;
- };
-
- 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);
- return res;
-}
-
// clang-format off
-DECLARE_DETELEMENT(ce_MRICH, createDetector)
+DECLARE_DETELEMENT(refdet_ce_MRICH, createDetector)
diff --git a/src/ref_utils.cpp b/src/ref_utils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..629679bdc2acdbe217e3322b45d08225b973dcfa
--- /dev/null
+++ b/src/ref_utils.cpp
@@ -0,0 +1,86 @@
+#include "ref_utils.h"
+
+// some utility functions that can be shared
+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)
+{
+ if (pt.r() > rmax || pt.r() < rmin) {
+ return false;
+ }
+
+ // 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.),
+ };
+ for (auto &p : pts) {
+ if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
+ return false;
+ }
+ }
+ return true;
+}
+
+// check if a square is overlapped with the others
+inline bool overlap(const Point &pt, double side, 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)) {
+ return true;
+ }
+ }
+ return false;
+}
+
+// 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)
+{
+ // outside of the ring or overlapping
+ if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, 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);
+}
+
+// fill squares
+std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax, double phmin, double phmax)
+{
+ // 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);
+
+ auto find_seed = [] (const Point &ref, int n, double side, 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)) {
+ return pt;
+ }
+ }
+ }
+ return ref;
+ };
+
+ 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);
+ return res;
+}
+
+
+} // ref::utils
diff --git a/src/ref_utils.h b/src/ref_utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac14e4df2ca32b2d9009544d0cedef6deea314de
--- /dev/null
+++ b/src/ref_utils.h
@@ -0,0 +1,14 @@
+#pragma once
+#include <vector>
+#include "Math/Point2D.h"
+
+// some utility functions that can be shared
+namespace ref::utils {
+
+typedef ROOT::Math::XYPoint Point;
+
+// fill squares in a ring
+std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax,
+ double phmin = 0., double phmax = 2.*M_PI);
+
+} // ref::utils