Commit a1fb656e authored by Chao Peng's avatar Chao Peng
Browse files

Resolve "Implement ce_MRICH"

parent 7772361a
<?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>
......@@ -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>
......
......@@ -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"/>
-->
......
......@@ -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.));
......
#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)
#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
#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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment