Skip to content
Snippets Groups Projects

Detailed Forward RICH Geometry

Merged Detailed Forward RICH Geometry
(removed):master into master
Merged Chao Peng requested to merge (removed):master into master
Compare
and
3 files
+ 275
− 18
Compare changes
  • Side-by-side
  • Inline
Files
3
src/ForwardRICH_geo.cpp 0 → 100644
+ 236
− 0
//==========================================================================
// 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)
{
// 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
double halfLength = 0.5 * (mirZ + 5.0*cm);
double rmin = std::min(pRmin, tRmin);
double rmax = std::max(std::max(pRmax, tRmax2), tRmax1);
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; }
}
Tube envShape(std::max(0., rmin - 0.1*cm), rmax + 0.1*cm, halfLength, 0., 2*M_PI);
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 radius = sl.attr<double>(_Unicode(curve));
auto mRmin = sl.attr<double>(_Unicode(rmin));
auto mRmax = sl.attr<double>(_Unicode(rmax));
double mTheta1 = std::asin(mRmin/radius);
double mTheta2 = std::asin(mRmax/radius);
double rotY = -std::asin(focus/radius);
// mirror slice shape
// somehow geant4 does not support -wphi/2. to wphi/2., so additonal rotation in Z
Sphere mirShape(radius, radius + mThick, mTheta1, mTheta2, 0., wphi);
Volume mirVol(Form("mirror_v_dummy%d", imod), mirShape, mirMat);
mirVol.setVisAttributes(desc.visAttributes(mir.visStr()));
// action is in a reverse order
Transform3D tr = Translation3D(0., 0., mirZ - halfLength)
* RotationZ(rotZ)
* RotationY(rotY)
* Translation3D(0., 0., -radius)
* RotationZ(-wphi/2.);
DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
auto mirPV = envVol.placeVolume(mirVol, tr);
mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", 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)