Skip to content
Snippets Groups Projects
Select Git revision
  • beampipe_overlap
  • master default protected
  • 24-optimize-barrel-em-calo-geometry
  • 18-make-ecal-blue
  • new_runner
  • geo_test
  • view_test
  • bigger_solenoid
  • vertex_tracker_radius
  • move_rich
  • rp_overlap
  • rich_optimize
  • acts_v4
  • tracker2
  • rich_test
  • no_rich_hcal
  • dawnfix
  • ecal_radius
  • v0.0.2
  • v0.0.1
20 results

ForwardRICH_geo.cpp

Blame
  • Chao Peng's avatar
    Chao Peng authored
    a6b04249
    History
    ForwardRICH_geo.cpp 9.76 KiB
    //==========================================================================
    //  Forward Ring Imaging Cherenkov Detector
    //--------------------------------------------------------------------------
    //
    // Author: C. Peng (ANL)
    // Date: 09/30/2020
    //
    //==========================================================================
    
    #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"
    
    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)
    {
        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) {
                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)
    {
        // outside of the ring or overlapping
        if (!in_ring(p, lside, rmin, rmax) || 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);
    }
    
    // fill squares
    std::vector<Point> fill_squares(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);
    
        auto find_seed = [] (const Point &ref, int n, double side, double rmin, double rmax) {
            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)) {
                        return pt;
                    }
                }
            }
            return ref;
        };
    
        std::vector<Point> res;
        ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax);
        add_square(ref, res, lside, rmin, rmax);
        return res;
    }
    
    // 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 tank = detElem.child(_Unicode(tank));
        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));
    
        // mirror setting
        auto mThick = mir.thickness();
        auto mirZ = mir.attr<double>(_Unicode(zdiff));
    
        // 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(zdiff));
    
        // materials
        auto mirMat = desc.material(mir.materialStr());
        auto gasMat = desc.material(tank.attr<std::string>(_Unicode(gas)));
        auto mcpMat = desc.material(mcp.materialStr());
    
        // 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);
    
        UnionSolid envShape(env1, env2, Position(0., 0., -halfLength + pZ + pThick/2.));
    
        Volume envVol(detName + "_envelope", envShape, desc.material("AirOptical"));
        // envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
    
        // ---------------
        // Gas radiator container and 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();
        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, mTheta2, 0., wphi));
                // action is in a reverse order
                Transform3D tr = Translation3D(0., 0., mirZ - halfLength)   // move for z position
                               * RotationZ(rotZ)                            // rotate phi angle
                               * RotationY(rotY)                            // rotate for focus point
                               * 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 - 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.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++;
    
        // ---------------
        // 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);
    
        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);
        }
    
        Volume motherVol = desc.pickMotherVolume(det);
        PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0 + halfLength));
        envPV.addPhysVolID("system", detID);
        det.setPlacement(envPV);
    
        return det;
    }
    //@}
    
    // clang-format off
    DECLARE_DETELEMENT(ForwardRICH, createDetector)