Commit 5363fb42 authored by Chao Peng's avatar Chao Peng Committed by Whitney Armstrong
Browse files

move detector planes

Merge gas tank into the envelope

adjust envelope shape

Fixed the ID for the readout
parent 30fe162b
......@@ -25,7 +25,7 @@ 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)
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;
......@@ -39,7 +39,7 @@ inline bool in_ring(const Point &pt, double side, double rmin, double rmax)
Point(pt.x() + side/2., pt.y() + side/2.),
};
for (auto &p : pts) {
if (p.r() > rmax || p.r() < rmin) {
if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
return false;
}
}
......@@ -59,34 +59,36 @@ inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts)
}
// 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_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) || overlap(p, lside, res)) {
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);
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);
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)
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) {
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)) {
if (in_ring(pt, side, rmin, rmax, phmin, phmax)) {
return pt;
}
}
......@@ -95,8 +97,8 @@ std::vector<Point> fill_squares(Point ref, double lside, double rmin, double rma
};
std::vector<Point> res;
ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax);
add_square(ref, res, lside, rmin, rmax);
ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax, phmin, phmax);
add_square(ref, res, lside, rmin, rmax, phmin, phmax);
return res;
}
......@@ -110,19 +112,16 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
DetElement det(detName, detID);
xml::Component dims = detElem.dimensions();
xml::Component tank = detElem.child(_Unicode(tank));
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();
// 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));
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();
......@@ -139,42 +138,30 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
// materials
auto mirMat = desc.material(mir.materialStr());
auto gasMat = desc.material(tank.attr<std::string>(_Unicode(gas)));
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 = 0.5 * (mirZ + 5.0*cm);
double rmin = tRmin;
double rmax = tRmax2;
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; }
}
rmin = std::max(0., rmin - 0.1*cm);
rmax = 2.*halfLength/mirZ*rmax;
Cone env1(halfLength, rmin, tRmax1 + 0.1*cm , rmin, rmax + 0.1*cm);
// disk for detection plane
Tube env2(std::max(0., pRmin - 0.1*cm), pRmax + pTol + pGap + 0.1*cm, pThick + 0.1*cm, 0., 2.*M_PI);
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., -halfLength + pZ + pThick/2.));
UnionSolid envShape(env1, env2, Position(0., 0., pZ));
Volume envVol(detName + "_envelope", envShape, desc.material("AirOptical"));
// envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
Volume envVol(detName + "_envelope", envShape, gasMat);
envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
// ---------------
// Gas radiator container and spherical mirrors inside it
// 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();
......@@ -207,6 +194,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
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);
......@@ -216,7 +204,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
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 = envVol.placeVolume(mirVol, tr);
}
mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", imod);
DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
......@@ -227,6 +215,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
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);
......@@ -240,15 +229,26 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
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);
// 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);
......
......@@ -5,18 +5,19 @@
<detectors>
<detector id="ForwardRICH_ID" name="ForwardRICH" type="ForwardRICH" readout="ForwardRICHHits" vis="BlueVis">
<dimensions z0="RICHZMin" />
<mcppmt zdiff="1.0*cm" rmin="RICHRMin" rmax="RICHRMin+60*cm" rtol="1.0*cm" 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="50*cm" curve="0" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="0*degree" />
<slice focus="50*cm" curve="0" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="60*degree" />
<slice focus="50*cm" curve="0" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="120*degree" />
<slice focus="50*cm" curve="0" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="180*degree" />
<slice focus="50*cm" curve="0" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="240*degree" />
<slice focus="50*cm" curve="0" rmin="RICHRMin" rmax="RICHRMin+80*cm" phiw="59*degree" rotz="300*degree" />
<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>
......@@ -24,7 +25,7 @@
<readouts>
<readout name="ForwardRICHHits">
<segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" />
<id>system:5,layer:4,module:14,x:32:-16,y:-16</id>
<id>system:8,layer:4,piece:4,module:14,x:32:-16,y:-16</id>
</readout>
</readouts>
......
Markdown is supported
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