Improve RICH Design

The goal is to find a design that fits the need of EIC. It may include:

1. change of depth/gas of the tank
2. change of position/angle/size of the detection planes
3. change of position/angle/size/curvature of the mirrors

src/ForwardRICH_geo.cpp

@@ -25,7 +25,7 @@ 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)
+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;
@@ -39,7 +39,7 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax)
         Point(pt.x() + side/2., pt.y() + side/2.),
     };
     for (auto &p : pts) {
-        if (p.r() > rmax || p.r() < rmin) {
+        if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
             return false;
         }
     }
@@ -59,34 +59,36 @@ 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)
+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) || overlap(p, lside, res)) {
+    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);
-    add_square(Point(p.x() - lside, p.y()), res, lside, rmin, rmax);
-    add_square(Point(p.x(), p.y() + lside), res, lside, rmin, rmax);
-    add_square(Point(p.x(), p.y() - lside), res, lside, rmin, rmax);
+    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)
+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) {
+    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)) {
+                if (in_ring(pt, side, rmin, rmax, phmin, phmax)) {
                     return pt;
                 }
             }
@@ -95,8 +97,8 @@ std::vector<Point> fill_squares(Point ref, double lside, double rmin, double rma
     };
 
     std::vector<Point> res;
-    ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax);
-    add_square(ref, res, lside, rmin, rmax);
+    ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax, phmin, phmax);
+    add_square(ref, res, lside, rmin, rmax, phmin, phmax);
     return res;
 }
 
@@ -110,19 +112,16 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
 
     DetElement det(detName, detID);
     xml::Component dims = detElem.dimensions();
-    xml::Component tank = detElem.child(_Unicode(tank));
+    xml::Component rads = detElem.child(_Unicode(radiator));
     xml::Component mir = detElem.child(_Unicode(mirror));
     xml::Component mcp = detElem.child(_Unicode(mcppmt));
 
     // dimensions
     double z0 = dims.z0();
-
-    // gas tank
-    auto tRmin = tank.rmin();
-    auto tRmax1 = tank.rmax1();
-    auto tRmax2 = tank.rmax2();
-    auto tLength = tank.length();
-    auto tZ = tank.attr<double>(_Unicode(zdiff));
+    double length = dims.length();
+    double rmin = dims.rmin();
+    double rmax1 = dims.attr<double>(_Unicode(rmax1));
+    double rmax2 = dims.attr<double>(_Unicode(rmax2));
 
     // mirror setting
     auto mThick = mir.thickness();
@@ -139,42 +138,30 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
 
     // materials
     auto mirMat = desc.material(mir.materialStr());
-    auto gasMat = desc.material(tank.attr<std::string>(_Unicode(gas)));
+    auto gasMat = desc.material(rads.materialStr());
     auto mcpMat = desc.material(mcp.materialStr());
 
+    // constants
+    auto richCenterAngle = std::atan((rmin + (rmax2 - rmin)/2.)/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 = 0.5 * (mirZ + 5.0*cm);
-    double rmin = tRmin;
-    double rmax = tRmax2;
-    for (xml::Collection_t sl(mir, _Unicode(slice)); sl; ++sl) {
-        auto mRmin = sl.attr<double>(_Unicode(rmin));
-        auto mRmax = sl.attr<double>(_Unicode(rmax));
-        if (mRmin < rmin) { rmin = mRmin; }
-        if (mRmax > rmax) { rmax = mRmax; }
-    }
-    rmin = std::max(0., rmin - 0.1*cm);
-    rmax = 2.*halfLength/mirZ*rmax;
-    Cone env1(halfLength, rmin, tRmax1 + 0.1*cm , rmin, rmax + 0.1*cm);
-    // disk for detection plane
-    Tube env2(std::max(0., pRmin - 0.1*cm), pRmax + pTol + pGap + 0.1*cm, pThick + 0.1*cm, 0., 2.*M_PI);
+    double halfLength = length/2.;
+    Cone env1(halfLength, rmin, rmax1, rmin, rmax2);
+    // envelope for detection plane
+    // Cone env2(halfLength - pZ/2., rmin, pRmax, rmin, rmax2);
+    Tube env2(rmin, pRmax + pTol + pGap + 1.0*cm, (length - pZ)/2., 0., 2*M_PI);
 
-    UnionSolid envShape(env1, env2, Position(0., 0., -halfLength + pZ + pThick/2.));
+    UnionSolid envShape(env1, env2, Position(0., 0., pZ));
 
-    Volume envVol(detName + "_envelope", envShape, desc.material("AirOptical"));
-    // envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
+    Volume envVol(detName + "_envelope", envShape, gasMat);
+    envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
 
     // ---------------
-    // Gas radiator container and spherical mirrors inside it
+    // spherical mirrors inside it
     int ilayer = 1;
-    Cone tankShape(tLength/2.0, tRmin, tRmax1, tRmin, tRmax2);
-    Volume tankVol("RICH_tank", tankShape, gasMat);
-    tankVol.setVisAttributes(desc.visAttributes(tank.visStr()));
-    auto tankPV = envVol.placeVolume(tankVol, Position(0., 0., -halfLength + tZ + tLength/2.0));
-    tankPV.addPhysVolID("layer", ilayer++);
-    DetElement tankDE(det, "Tank_DE", 1);
-    tankDE.setPlacement(tankPV);
 
     // optical surface
     OpticalSurfaceManager surfMgr = desc.surfaceManager();
@@ -207,6 +194,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
             Transform3D tr = Translation3D(0., 0., mirZ - halfLength)   // 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);
@@ -216,7 +204,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
             Transform3D tr = Translation3D(0., 0., mirZ - halfLength)   // 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 = envVol.placeVolume(mirVol, tr);
         }
         mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", imod);
         DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
@@ -227,6 +215,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
     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);
@@ -240,15 +229,26 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
     sens.setType("photoncounter");
     mcpVol.setSensitiveDetector(sens);
 
-    auto points = fill_squares(Point(0., 0.), pSize + pGap, pRmin - pTol - pGap, pRmax + pTol + pGap);
-    for (size_t i = 0; i < points.size(); ++i) {
-        auto pt = points[i];
-        auto mcpPV = envVol.placeVolume(mcpVol, Position(pt.x(), pt.y(), -halfLength + pZ + pThick/2.0));
-        mcpPV.addPhysVolID("layer", ilayer).addPhysVolID("module", i + 1);
-        DetElement mcpDE(det, Form("MCPPMT_DE%d", i + 1), i + 1);
-        mcpDE.setPlacement(mcpPV);
+    // photo-detector plane envelope
+    for (size_t ipd = 0; ipd < 6; ++ipd) {
+        double phmin = -M_PI/6.;
+        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);
+        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., -halfLength + 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 + halfLength));
     envPV.addPhysVolID("system", detID);