Skip to content
Snippets Groups Projects
Commit e7a63ba1 authored by Chao Peng's avatar Chao Peng
Browse files

Geometry update for Forward RICH detector

parent 327a35c0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
athena.xml
+ 1
− 1
View file @ e7a63ba1
...@@ -190,7 +190,7 @@ ...@@ -190,7 +190,7 @@
<include ref="compact/ce_mrich.xml"/> <include ref="compact/ce_mrich.xml"/>
<include ref="compact/tof_endcap.xml"/> <include ref="compact/tof_endcap.xml"/>
<include ref="compact/forward_trd.xml"/> <include ref="compact/forward_trd.xml"/>
<include ref="compact/forward_rich.xml"/> <include ref="compact/gaseous_rich.xml"/>
<include ref="ip6/B0_tracker.xml"/> <include ref="ip6/B0_tracker.xml"/>
<include ref="ip6/far_forward_offM_tracker.xml"/> <include ref="ip6/far_forward_offM_tracker.xml"/>
......
compact/forward_rich.xml deleted 100644 → 0
+ 0
− 68
View file @ 327a35c0
<?xml version="1.0" encoding="UTF-8"?>
<lccdd>
<define>
<constant name="ForwardRICH_zmin" value="BarrelTracking_length/2.0 + ForwardTracking_length "/>
<constant name="ForwardRICH_rmin" value="ForwardPID_rmin1"/>
<constant name="ForwardRICH_rmax0" value="TrackerBarrel_rmax"/>
<constant name="ForwardRICH_rmax1" value="EcalBarrel_rmin"/>
<constant name="ForwardRICH_rmax2" value="Solenoid_rmin-6*cm"/>
<constant name="ForwardRICHDepth" value="0.9*m"/>
</define>
<detectors>
<detector
id="ForwardRICH_ID"
name="ForwardRICH"
type="refdet_ForwardRICH"
readout="ForwardRICHHits"
vis="RICHVis">
<dimensions
z0="ForwardRICH_zmin"
snout_length="ForwardRICH_length - ForwardRICHDepth"
length="ForwardRICH_length"
rmin="ForwardRICH_rmin"
rmax0="ForwardRICH_rmax0"
rmax1="ForwardRICH_rmax1"
rmax2="ForwardRICH_rmax2"/>
<radiator material="N2cherenkov" />
<tank
length="ForwardRICHDepth"
gas="N2cherenkov"
vis="GreenVis"
rmin="ForwardRICH_rmin"
rmax1="ForwardRICH_rmax1"
rmax2="ForwardRICH_rmax2" />
<comment> What are the following MCP-PMT parameters?</comment>
<mcppmt
z0="0.0*cm"
rmin="ForwardRICH_rmax1"
rmax="ForwardRICH_rmax2"
rtol="1.0*cm"
vis="BlueVis"
module_size="10*cm"
module_gap="0.2*cm"
thickness="1.0*cm"
material="Quartz" />
<mirror
z0="ForwardRICH_length - 40*cm"
thickness="1*mm"
material="PyrexGlass"
vis="GrayVis">
<slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="0*degree" />
<slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="60*degree" />
<slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="120*degree" />
<slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="180*degree" />
<slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" phiw="59*degree" rotz="240*degree" />
<slice focus="10*cm" curve="300*cm" rmin="ForwardRICH_rmin" rmax="ForwardRICH_rmax1" 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>
compact/gaseous_rich.xml 0 → 100644
+ 73
− 0
View file @ e7a63ba1
<?xml version="1.0" encoding="UTF-8"?>
<lccdd>
<define>
<constant name="ForwardRICH_zmin" value="BarrelTracking_length/2.0 + ForwardTracking_length "/>
<constant name="ForwardRICH_rmin" value="ForwardPID_rmin1"/>
<constant name="ForwardRICH_rmax0" value="TrackerBarrel_rmax"/>
<constant name="ForwardRICH_rmax1" value="EcalBarrel_rmin"/>
<constant name="ForwardRICH_rmax2" value="Solenoid_rmin-6*cm"/>
<constant name="ForwardRICHDepth" value="0.9*m"/>
<constant name="ForwardRICH_sensor_rmin" value="ForwardRICH_rmax1 - 40.*cm"/>
<constant name="ForwardRICH_sensor_rmax" value="ForwardRICH_rmax2"/>
</define>
<detectors>
<detector
id="ForwardRICH_ID"
name="GaseousRICH"
type="athena_GaseousRICH"
readout="ForwardRICHHits"
gas="N2cherenkov"
vis="RICHVis">
<dimensions
z0="ForwardRICH_zmin"
snout_length="ForwardRICH_length - ForwardRICHDepth"
length="ForwardRICH_length"
rmin="ForwardRICH_rmin"
rmax0="ForwardRICH_rmax0"
rmax1="ForwardRICH_rmax1"
rmax2="ForwardRICH_rmax2"/>
<mirrors thickness="1*mm" material="PyrexGlass" vis="GrayVis">
<position z="ForwardRICH_length - 13*cm"/>
<mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
rotz="0*degree" roty="-15*degree" rotx="0"/>
<mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
rotz="60*degree" roty="-15*degree" rotx="0"/>
<mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
rotz="120*degree" roty="-15*degree" rotx="0"/>
<mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
rotz="180*degree" roty="-15*degree" rotx="0"/>
<mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
rotz="240*degree" roty="-15*degree" rotx="0"/>
<mirror curve="300*cm" rmin="ForwardRICH_rmin+5*mm" rmax="ForwardRICH_rmax1" phiw="58*degree"
rotz="300*degree" roty="-15*degree" rotx="0"/>
</mirrors>
<sensors>
<module sx="10*cm" sy="10*cm" sz="1.0*cm" gap="0.2*cm" material="Quartz" vis="AnlGold"/>
<comment>A thin optical material to accept optical photon in simulation</comment>
<optical material="AirOptical" thickness="0.1*mm"/>
<position z="8.*cm"/>
<sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
rotz="0*degree" roty="-15*degree" rotx="0"/>
<sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
rotz="60*degree" roty="-15*degree" rotx="0"/>
<sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
rotz="120*degree" roty="-15*degree" rotx="0"/>
<sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
rotz="180*degree" roty="-15*degree" rotx="0"/>
<sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
rotz="240*degree" roty="-15*degree" rotx="0"/>
<sector rmin="ForwardRICH_sensor_rmin" rmax="ForwardRICH_sensor_rmax" phiw="58*degree"
rotz="300*degree" roty="-15*degree" rotx="0"/>
</sensors>
</detector>
</detectors>
<readouts>
<readout name="ForwardRICHHits">
<segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" />
<id>system:8,sector:8,module:16,x:32:-16,y:-16</id>
</readout>
</readouts>
</lccdd>
src/ForwardRICH_geo.cpp deleted 100644 → 0
+ 0
− 192
View file @ 327a35c0
//==========================================================================
// Forward Ring Imaging Cherenkov Detector
//--------------------------------------------------------------------------
//
// Author: C. Peng (ANL)
// Date: 09/30/2020
//
//==========================================================================
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "GeometryHelpers.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;
// 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));
xml::Component tank = detElem.child(_Unicode(tank));
// dimensions
double z0 = dims.z0();
double length = dims.length();
double rmin = dims.rmin();
double rmax0 = dims.attr<double>(_Unicode(rmax0));
double rmax1 = dims.attr<double>(_Unicode(rmax1));
double rmax2 = dims.attr<double>(_Unicode(rmax2));
double snout_length = dims.attr<double>(_Unicode(snout_length));
// mirror setting
auto mThick = mir.thickness();
auto mirZ = mir.attr<double>(_Unicode(z0));
// 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(z0));
// tank parameters
double tank_length = length - snout_length;
// materials
auto mirMat = desc.material(mir.materialStr());
auto gasMat = desc.material(rads.materialStr());
auto mcpMat = desc.material(mcp.materialStr());
double front_offset = snout_length+tank_length/2.0;
// constants
auto richCenterAngle = std::atan((rmin + (rmax1 - rmin)/2.)/(front_offset+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(snout_length/2.0, rmin, rmax0, rmin, rmax1);
// envelope for detection plane
// Cone env2(halfLength - pZ/2., rmin, pRmax, rmin, rmax2);
Tube env2(rmin, pRmax + pTol + pGap + 1.0*cm, tank_length/2., 0., 2*M_PI);
UnionSolid envShape(env2, env1, Position(0., 0., -tank_length/2.-snout_length/2));
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*1.01, mTheta2*0.99, 0., wphi));
// action is in a reverse order
Transform3D tr = Translation3D(0., 0., mirZ - front_offset) // 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 - front_offset) // 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.5; // added 0.5 to make it smaller
double phmax = M_PI/6.5;
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 = 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.));
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., -front_offset + 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 + front_offset));
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
//@}
// clang-format off
DECLARE_DETELEMENT(refdet_ForwardRICH, createDetector)
src/GaseousRICH_geo.cpp 0 → 100644
+ 215
− 0
View file @ e7a63ba1
//==========================================================================
// Gaseous Ring Imaging Cherenkov Detector
//--------------------------------------------------------------------------
//
// Author: C. Peng (ANL)
// Date: 09/30/2020
//
//==========================================================================
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "GeometryHelpers.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;
// headers
void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset);
void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset);
void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens);
// helper function to get x, y, z if defined in a xml component
template<class XmlComp>
Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
{
Position pos(0., 0., 0.);
if (comp.hasChild(name)) {
auto child = comp.child(name);
pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
}
return pos;
}
// 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();
// build a big envelope for all the components, filled with optical material to ensure transport of optical photons
double z0 = dims.z0();
double length = dims.length();
double rmin = dims.rmin();
double rmax0 = dims.attr<double>(_Unicode(rmax0));
double rmax1 = dims.attr<double>(_Unicode(rmax1));
double rmax2 = dims.attr<double>(_Unicode(rmax2));
double snout_length = dims.attr<double>(_Unicode(snout_length));
// fill envelope with radiator materials (need optical property for optical photons)
auto gasMat = desc.material(dd4hep::getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
Cone snout(snout_length/2.0, rmin, rmax0, rmin, rmax1);
Tube tank(rmin, rmax2, (length - snout_length)/2., 0., 2*M_PI);
// shift the snout to the left side of tank
UnionSolid envShape(tank, snout, Position(0., 0., -length/2.));
// some reference points
auto snout_front = Position(0., 0., -(length + snout_length)/2.);
// auto snout_end = Position(0., 0., -(length - snout_length)/2.); // tank_front
// auto tank_end = Position(0., 0., (length - snout_length)/2.);
Volume envVol(detName + "_envelope", envShape, gasMat);
envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
// sensitive detector type
sens.setType("photoncounter");
// @TODO: place a radiator
// build_radiator(desc, envVol, detElement.child(_Unicode(radiator)), snout_front);
// place mirrors
build_mirrors(desc, det, envVol, detElem.child(_Unicode(mirrors)), snout_front);
// place photo-sensors
build_sensors(desc, envVol, detElem.child(_Unicode(sensors)), snout_front, sens);
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0) - snout_front);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// @TODO: implement a radiator, now the envelope serves as the radiator
void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset)
{
// place holder
}
// place mirrors
void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset)
{
double thickness = dd4hep::getAttrOrDefault<double>(plm, _Unicode(dz), 1.*dd4hep::mm);
auto mat = desc.material(plm.attr<std::string>(_Unicode(material)));
auto vis = desc.visAttributes(plm.attr<std::string>(_Unicode(vis)));
// optical surface
OpticalSurfaceManager surfMgr = desc.surfaceManager();
auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(plm, _Unicode(surface), "MirrorOpticalSurface"));
// placements
auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
auto grot = get_xml_xyz(plm, _Unicode(position));
int imod = 1;
for (xml::Collection_t mir(plm, _Unicode(mirror)); mir; ++mir, ++imod) {
double rmin = mir.attr<double>(_Unicode(rmin));
double rmax = mir.attr<double>(_Unicode(rmax));
double phiw = dd4hep::getAttrOrDefault<double>(mir, _Unicode(phiw), 2.*M_PI);
double rotz = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotz), 0.);
double roty = dd4hep::getAttrOrDefault<double>(mir, _Unicode(roty), 0.);
double rotx = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotx), 0.);
Volume vol(Form("mirror_v_dummy%d", imod));
vol.setMaterial(mat);
vol.setVisAttributes(vis);
// mirror curvature
double curve = dd4hep::getAttrOrDefault<double>(mir, _Unicode(curve), 0.);
// spherical mirror
if (curve > 0.) {
double th1 = std::asin(rmin/curve);
double th2 = std::asin(rmax/curve);
vol.setSolid(Sphere(curve, curve + thickness, th1, th2, 0., phiw));
// plane mirror
} else {
vol.setSolid(Tube(rmin, rmax, thickness/2., 0., phiw));
}
// transforms are in a reverse order
Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
* RotationZYX(grot.z(), grot.y(), grot.x())
* RotationZ(rotz) // rotation of the piece
* RotationY(roty) // rotation of the piece
* RotationX(rotx) // rotation of the piece
* Translation3D(0., 0., -curve) // move spherical shell to origin (no move for planes)
* RotationZ(-phiw/2.); // center phi angle to 0. (-phiw/2., phiw/2.)
auto pv = env.placeVolume(vol, tr);
DetElement de(sdet, Form("mirror_de%d", imod), imod);
de.setPlacement(pv);
SkinSurface skin(desc, de, Form("mirror_optical_surface%d", imod), surf, vol);
skin.isValid();
}
}
// place photo-sensors
void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens)
{
// build sensor unit geometry
auto mod = plm.child(_Unicode(module));
double sx = mod.attr<double>(_Unicode(sx));
double sy = mod.attr<double>(_Unicode(sy));
double sz = mod.attr<double>(_Unicode(sz));
double gap = mod.attr<double>(_Unicode(gap));
auto mat = desc.material(mod.attr<std::string>(_Unicode(material)));
auto vis = desc.visAttributes(mod.attr<std::string>(_Unicode(vis)));
Box sensor(sx/2., sy/2., sz/2.);
Volume svol("sensor_v", sensor, mat);
svol.setVisAttributes(vis);
// a thin layer of cherenkov gas for accepting optical photons
auto opt = plm.child(_Unicode(optical));
double opthick = opt.attr<double>(_Unicode(thickness));
auto opmat = desc.material(opt.attr<std::string>(_Unicode(material)));
Box opshape(sx/2., sy/2., sz/2. + opthick/2.);
Volume opvol("sensor_v_optical", opshape, opmat);
opvol.placeVolume(svol, Position(0., 0., 0.));
opvol.setSensitiveDetector(sens);
// photo-detector plane envelope
auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
auto grot = get_xml_xyz(plm, _Unicode(position));
int isec = 1;
for (xml::Collection_t sec(plm, _Unicode(sector)); sec; ++sec, ++isec) {
double rmin = sec.attr<double>(_Unicode(rmin));
double rmax = sec.attr<double>(_Unicode(rmax));
double phiw = dd4hep::getAttrOrDefault<double>(sec, _Unicode(phiw), 2.*M_PI);
double rotz = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotz), 0.);
double roty = dd4hep::getAttrOrDefault<double>(sec, _Unicode(roty), 0.);
double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
// fill sensors to the piece
auto points = ref::utils::fillRectangles({0., 0.}, sx + gap, sy + gap, rmin - gap, rmax + gap, -phiw/2., phiw/2.);
int imod = 1;
for (auto &p : points) {
// transofrms are in a reversed order
Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
* RotationZYX(grot.z(), grot.y(), grot.x())
* RotationZ(rotz) // rotation of the sector
* RotationY(roty) // rotation of the sector
* RotationX(rotx) // rotation of the sector
* Translation3D(p.x(), p.y(), 0.); // place modules in each sector
auto pv = env.placeVolume(opvol, tr);
pv.addPhysVolID("sector", isec).addPhysVolID("module", imod++);
}
}
}
//@}
// clang-format off
DECLARE_DETELEMENT(athena_GaseousRICH, createDetector)
src/GeometryHelpers.cpp
+ 21
− 22
View file @ e7a63ba1
...@@ -6,7 +6,7 @@ namespace ref::utils { ...@@ -6,7 +6,7 @@ namespace ref::utils {
typedef ROOT::Math::XYPoint Point; typedef ROOT::Math::XYPoint Point;
// check if a square in a ring // check if a square in a ring
inline bool in_ring(const Point &pt, double side, double rmin, double rmax, double phmin, double phmax) inline bool in_ring(const Point &pt, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{ {
if (pt.r() > rmax || pt.r() < rmin) { if (pt.r() > rmax || pt.r() < rmin) {
return false; return false;
...@@ -14,10 +14,10 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax, doub ...@@ -14,10 +14,10 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax, doub
// check four corners // check four corners
std::vector<Point> pts { std::vector<Point> pts {
Point(pt.x() - side/2., pt.y() - side/2.), Point(pt.x() - sx/2., pt.y() - sy/2.),
Point(pt.x() - side/2., pt.y() + side/2.), Point(pt.x() - sx/2., pt.y() + sy/2.),
Point(pt.x() + side/2., pt.y() - side/2.), Point(pt.x() + sx/2., pt.y() - sy/2.),
Point(pt.x() + side/2., pt.y() + side/2.), Point(pt.x() + sx/2., pt.y() + sy/2.),
}; };
for (auto &p : pts) { for (auto &p : pts) {
if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) { if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
...@@ -28,11 +28,11 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax, doub ...@@ -28,11 +28,11 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax, doub
} }
// check if a square is overlapped with the others // check if a square is overlapped with the others
inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts) inline bool overlap(const Point &pt, double sx, double sy, const std::vector<Point> &pts)
{ {
for (auto &p : pts) { for (auto &p : pts) {
auto pn = (p - pt)/side; auto pn = p - pt;
if ((std::abs(pn.x()) < 1. - 1e-6) && (std::abs(pn.y()) < 1. - 1e-6)) { if ((std::abs(pn.x()) < (1. - 1e-6)*sx) && (std::abs(pn.y()) < (1. - 1e-6)*sy)) {
return true; return true;
} }
} }
...@@ -40,25 +40,24 @@ inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts) ...@@ -40,25 +40,24 @@ inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts)
} }
// a helper function to recursively fill square in a ring // 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, void add_rectangle(Point p, std::vector<Point> &res, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
double phmin, double phmax)
{ {
// outside of the ring or overlapping // outside of the ring or overlapping
if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, res)) { if (!in_ring(p, sx, sy, rmin, rmax, phmin, phmax) || overlap(p, sx, sy, res)) {
return; return;
} }
res.emplace_back(p); res.emplace_back(p);
// check adjacent squares // check adjacent squares
add_square(Point(p.x() + lside, p.y()), res, lside, rmin, rmax, phmin, phmax); add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax);
add_square(Point(p.x() - lside, p.y()), res, lside, rmin, rmax, phmin, phmax); add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax);
add_square(Point(p.x(), p.y() + lside), res, lside, rmin, rmax, phmin, phmax); add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax);
add_square(Point(p.x(), p.y() - lside), res, lside, rmin, rmax, phmin, phmax); add_rectangle(Point(p.x(), p.y() - sy), res, sx, sy, rmin, rmax, phmin, phmax);
} }
// fill squares // fill squares
std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax, double phmin, double phmax) std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{ {
// convert (0, 2pi) to (-pi, pi) // convert (0, 2pi) to (-pi, pi)
if (phmax > M_PI) { if (phmax > M_PI) {
...@@ -67,13 +66,13 @@ std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax ...@@ -67,13 +66,13 @@ std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax
} }
// start with a seed square and find one in the ring // start with a seed square and find one in the ring
// move to center // move to center
ref = ref - Point(int(ref.x()/lside)*lside, int(ref.y()/lside)*lside); ref = ref - Point(int(ref.x()/sx)*sx, int(ref.y()/sy)*sy);
auto find_seed = [] (const Point &ref, int n, double side, double rmin, double rmax, double phmin, double phmax) { auto find_seed = [] (const Point &ref, int n, double sx, double sy, double rmin, double rmax, double phmin, double phmax) {
for (int ix = -n; ix < n; ++ix) { for (int ix = -n; ix < n; ++ix) {
for (int iy = -n; iy < n; ++iy) { for (int iy = -n; iy < n; ++iy) {
Point pt(ref.x() + ix*side, ref.y() + iy*side); Point pt(ref.x() + ix*sx, ref.y() + iy*sy);
if (in_ring(pt, side, rmin, rmax, phmin, phmax)) { if (in_ring(pt, sx, sy, rmin, rmax, phmin, phmax)) {
return pt; return pt;
} }
} }
...@@ -82,8 +81,8 @@ std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax ...@@ -82,8 +81,8 @@ std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax
}; };
std::vector<Point> res; std::vector<Point> res;
ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax, phmin, phmax); ref = find_seed(ref, int(rmax/sx) + 2, sx, sy, rmin, rmax, phmin, phmax);
add_square(ref, res, lside, rmin, rmax, phmin, phmax); add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax);
return res; return res;
} }
......
src/GeometryHelpers.h
+ 8
− 2
View file @ e7a63ba1
...@@ -7,8 +7,14 @@ namespace ref::utils { ...@@ -7,8 +7,14 @@ namespace ref::utils {
typedef ROOT::Math::XYPoint Point; typedef ROOT::Math::XYPoint Point;
// fill rectangles in a ring
std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax,
double phmin = -M_PI, double phmax = M_PI);
// fill squares in a ring // fill squares in a ring
std::vector<Point> fillSquares(Point ref, double lside, double rmin, double rmax, inline std::vector<Point> fillSquares(Point ref, double size, double rmin, double rmax,
double phmin = -M_PI, double phmax = M_PI); double phmin = -M_PI, double phmax = M_PI)
{
return fillRectangles(ref, size, size, rmin, rmax, phmin, phmax);
}
} // ref::utils } // ref::utils
src/HomogeneousCalorimeter_geo.cpp
+ 1
− 1
View file @ e7a63ba1
...@@ -275,7 +275,7 @@ static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, Sens ...@@ -275,7 +275,7 @@ static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, Sens
double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.); double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI); double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
auto points = ref::utils::fillSquares({0., 0.}, modSize.x(), rmin, rmax, phimin, phimax); auto points = ref::utils::fillRectangles({0., 0.}, modSize.x(), modSize.y(), rmin, rmax, phimin, phimax);
// placement to mother // placement to mother
auto pos = get_xml_xyz(plm, _Unicode(position)); auto pos = get_xml_xyz(plm, _Unicode(position));
auto rot = get_xml_xyz(plm, _Unicode(rotation)); auto rot = get_xml_xyz(plm, _Unicode(rotation));
......
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