Merge branch '12-implement-ce_mrich' into 'master'

Resolve "Implement ce_MRICH"

Closes #12

See merge request !10

compact/ce_mrich.xml

+<?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>

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="-EcalEndcap_zmin+10.*cm"/>
+
+
 
 
   </define>

reference_detector.xml

@@ -119,6 +119,7 @@
   <include ref="compact/ecal.xml"/>
   <include ref="compact/hcal.xml"/>
   <include ref="compact/forward_rich.xml"/>
+  <include ref="compact/ce_mrich.xml"/>
   <!--
   <include ref="compact/roman_pots.xml"/>
   -->

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.));

src/ce_MRICH.cpp

+
+#include <XML/Helper.h>
+#include "TMath.h"
+#include "TString.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;
+
+
+void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, SensitiveDetector &sens);
+
+
+// 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));
+
+    auto rmin = dims.rmin();
+    auto rmax = dims.rmax();
+    auto length = dims.length();
+    auto z0 = dims.z();
+
+    auto gasMat = desc.material("AirOptical");
+
+    // detector envelope
+    Tube envShape(rmin, rmax, length/2., 0., 2*M_PI);
+    Volume envVol("ce_MRICH_GVol", envShape, gasMat);
+    envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
+
+    // modules
+    addModules(envVol, detElem, desc, sens);
+
+    // place envelope
+    Volume motherVol = desc.pickMotherVolume(det);
+    PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0));
+    envPV.addPhysVolID("system", detID);
+    det.setPlacement(envPV);
+    return det;
+}
+
+
+void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, SensitiveDetector &sens)
+{
+    xml::Component dims = detElem.dimensions();
+    xml::Component mods = detElem.child(_Unicode(modules));
+
+    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 modMat = desc.material(mods.materialStr());
+    auto gasMat = desc.material("AirOptical");
+
+    // single module
+    Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm);
+    Volume mVol("ce_MRICH_mod_Solid", mShape, modMat);
+
+    // a thin gas layer to detect optical photons
+    Box modShape(mWidth/2., mWidth/2., mThick/2.);
+    Volume modVol("ce_MRICH_mod_Solid_v", modShape, gasMat);
+    // thin gas layer is on top (+z) of the material
+    modVol.placeVolume(mVol, Position(0., 0., -0.1*mm));
+
+    modVol.setVisAttributes(desc.visAttributes(mods.visStr()));
+    sens.setType("photoncounter");
+    modVol.setSensitiveDetector(sens);
+
+    // place modules in the sectors (disk)
+    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.;
+    int imod = 1;
+    for (auto &p : points) {
+        // operations are inversely ordered
+        Transform3D tr = Translation3D(p.x(), p.y(), 0.)        // move to position
+                       * RotationY(roty);                       // facing z = 0.
+        auto modPV = mother.placeVolume(modVol, tr);
+        modPV.addPhysVolID("sector", 1).addPhysVolID("module", imod ++);
+    }
+}
+
+// clang-format off
+DECLARE_DETELEMENT(refdet_ce_MRICH, createDetector)
+

src/ref_utils.cpp

+#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

src/ref_utils.h

+#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 = -M_PI, double phmax = M_PI);
+
+} // ref::utils