Skip to content
Snippets Groups Projects
Commit 3803cbee authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Copied topside rich to reference_detector 

	new file:   compact/forward_rich.xml
	modified:   reference_detector.xml
	new file:   src/ForwardRICH_geo.cpp
parent a265ee2e
No related branches found
No related tags found
1 merge request!20Resolve "Implement Forward Rich"
Pipeline #7594 failed
Stage: config
Stage: build
Stage: docs
Stage: test
Stage: finalize
Stage: deploy
This commit is part of merge request !20. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
compact/forward_rich.xml 0 → 100644
+ 32
− 0
View file @ 3803cbee
<?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>
reference_detector.xml
+ 1
− 1
View file @ 3803cbee
......@@ -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"/>
-->
<!--
......
src/ForwardRICH_geo.cpp 0 → 100644
+ 263
− 0
View file @ 3803cbee
//==========================================================================
// 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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment