diff --git a/compact/forward_rich.xml b/compact/forward_rich.xml new file mode 100644 index 0000000000000000000000000000000000000000..e129e31442127c8c958719955f4c34f465785f95 --- /dev/null +++ b/compact/forward_rich.xml @@ -0,0 +1,32 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd> + <define> + </define> + + <detectors> + <detector id="ForwardRICH_ID" name="ForwardRICH" type="ForwardRICH" readout="ForwardRICHHits" vis="BlueVis"> + <dimensions z0="RICHZMin" length="RICHDepth+20*cm" rmin="RICHRMin" rmax1="RICHRMin+40*cm" rmax2="RICHRMin+80*cm"/> + <radiator material="N2cherenkov" /> + <mcppmt zdiff="15.0*cm" rmin="SolenoidYokeEndcapP_rmin - 10*cm" rmax="SolenoidYokeEndcapP_rmin + 80*cm" rtol="1.0*cm" vis="BlueVis" + module_size="10*cm" module_gap="0.2*cm" thickness="1.0*cm" material="Quartz" /> + <tank zdiff="5.0*cm" length="RICHDepth" gas="N2cherenkov" vis="GreenVis" + rmin="RICHRMin" rmax1="RICHRMin+40*cm" rmax2="RICHRMin+80*cm" /> + <mirror zdiff="RICHDepth+7.0*cm" thickness="1*mm" material="PyrexGlass" vis="GrayVis"> + <slice focus="10*cm" curve="300*cm" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="0*degree" /> + <slice focus="10*cm" curve="300*cm" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="60*degree" /> + <slice focus="10*cm" curve="300*cm" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="120*degree" /> + <slice focus="10*cm" curve="300*cm" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="180*degree" /> + <slice focus="10*cm" curve="300*cm" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="240*degree" /> + <slice focus="10*cm" curve="300*cm" rmin="RICHRMin" rmax="RICHRMin+80*cm" 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/reference_detector.xml b/reference_detector.xml index 2b2aa0a187a9972863582407d01c70a05726f43a..463faf881db2a2176313134253c3b407b009a4ca 100644 --- a/reference_detector.xml +++ b/reference_detector.xml @@ -117,8 +117,8 @@ <include ref="compact/solenoid.xml"/> <include ref="compact/ecal.xml"/> <include ref="compact/hcal.xml"/> + <include ref="compact/forward_rich.xml"/> <!-- - <include ref="compact/reference_detector_rich.xml"/> <include ref="compact/roman_pots.xml"/> --> <!-- diff --git a/src/ForwardRICH_geo.cpp b/src/ForwardRICH_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6520c30aac30492a8950c1590661bce57e089e47 --- /dev/null +++ b/src/ForwardRICH_geo.cpp @@ -0,0 +1,263 @@ +//========================================================================== +// 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, 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) +{ + 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)); + + // dimensions + double z0 = dims.z0(); + 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(); + 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(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 = 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., pZ)); + + 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, 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 + * 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 - 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++; + + // --------------- + // 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.; + 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); + det.setPlacement(envPV); + + return det; +} +//@} + +// clang-format off +DECLARE_DETELEMENT(refdet_ForwardRICH, createDetector) +