diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b0fff8df3b3dc0d8ab6e6286e2af5931f2e6ccfd..a95f75e2f6d461204cdcb2a1b0e50cbb1c3afca1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -12,7 +12,7 @@ default:
before_script:
- mkdir -p images && mkdir -p doc/
- git clone https://eicweb.phy.anl.gov/EIC/detectors/accelerator.git && ln -s accelerator/eic
- - git clone https://eicweb.phy.anl.gov/EIC/detectors/ip6.git && mkdir ip6_build && cd ip6_build && cmake ../ip6/. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 && make install && cd .. && ln -s ip6/ip6
+ - git clone https://eicweb.phy.anl.gov/EIC/detectors/ip6.git eic_ip6 && mkdir ip6_build && cd ip6_build && cmake ../eic_ip6/. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 && make install && cd .. && ln -s eic_ip6/ip6
- mkdir build && cd build && cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j20 && make install && cd ..
artifacts:
paths:
diff --git a/compact/definitions.xml b/compact/definitions.xml
index 750a8ab3309436f4690929947c059bdca708dae9..c7321b9729a1c651422d86b0627449cf6cac8872 100644
--- a/compact/definitions.xml
+++ b/compact/definitions.xml
@@ -423,6 +423,7 @@
</comment>
<constant name="RICHZMin" value="SolenoidYokeEndcap_zmin + 1 * cm"/>
<constant name="RICHRMin" value="15 * cm"/>
+ <constant name="RICHRMax" value="60 * cm"/>
<constant name="RICHDepth" value="1 * m"/>
diff --git a/compact/display.xml b/compact/display.xml
index b2ad9891af9b6c0d03f7b08abd2cd7becf420054..ae6af1bf18fd1819d37492fc0e8484eee122e111 100644
--- a/compact/display.xml
+++ b/compact/display.xml
@@ -53,11 +53,11 @@
<comment>
Deprecated colors.
</comment>
- <vis name="GreenVis" alpha="0.2" r= "0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="GreenVis" alpha="1.0" r= "0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
<vis name="RedVis" alpha="0.2" r= "1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/>
<vis name="RPVis" alpha="0.99" r= "1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/>
<vis name="RPLayerVis" alpha="0.99" r= "0.0" g="0.7" b="0.3" showDaughters="true" visible="true" lineStyle="solid" drawingStyle="solid" />
- <vis name="BlueVis" alpha="0.2" r= "0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
+ <vis name="BlueVis" alpha="1.0" r= "0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
<vis name="OrangeVis" alpha="1.0" r= "1.0" g="0.45" b="0.0" showDaughters="true" visible="true"/>
<vis name="RedGreenVis" alpha="0.5" r= "1.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
<vis name="BlueGreenVis" alpha="0.5" r= "0.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
@@ -66,5 +66,5 @@
<vis name="RBG015" alpha="0.5" r= "0.0" g=".2" b="1.0" showDaughters="true" visible="true"/>
<vis name="RBG510" alpha="0.5" r= "1.0" g=".2" b="0.0" showDaughters="true" visible="true"/>
<vis name="RBG" alpha="0.5" r= "1.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
- <vis name="GrayVis" alpha="0.5" r= "0.75" g="0.75" b="0.75" showDaughters="true" visible="true"/>
+ <vis name="GrayVis" alpha="1.0" r= "0.75" g="0.75" b="0.75" showDaughters="true" visible="true"/>
</display>
diff --git a/compact/forward_rich.xml b/compact/forward_rich.xml
new file mode 100644
index 0000000000000000000000000000000000000000..d94e9e162592d826ae9e7f1f14216cf6f39fc94f
--- /dev/null
+++ b/compact/forward_rich.xml
@@ -0,0 +1,32 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+ <define>
+ </define>
+
+ <detectors>
+ <detector id="ForwardRICH_ID" name="ForwardRICH" type="refdet_ForwardRICH" readout="ForwardRICHHits" vis="BlueVis">
+ <dimensions z0="RICHZMin" length="RICHDepth+20*cm" rmin="RICHRMin" rmax1="RICHRMax*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>
diff --git a/reference_detector.xml b/reference_detector.xml
index 2b2aa0a187a9972863582407d01c70a05726f43a..463faf881db2a2176313134253c3b407b009a4ca 100644
--- a/reference_detector.xml
+++ b/reference_detector.xml
@@ -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"/>
-->
<!--
diff --git a/scripts/view1/.DAWN_1.history b/scripts/view1/.DAWN_1.history
index 42729c65bc406cdb4b98cfef27a6e803548a0bba..3a43b62a9d1b3c8407d8d42151e4218e629c4f3c 100644
--- a/scripts/view1/.DAWN_1.history
+++ b/scripts/view1/.DAWN_1.history
@@ -5,10 +5,10 @@
0
0
0
-16.4
+1
1
0.001
-0
+1
1
1
1
diff --git a/scripts/view1/generate_eps b/scripts/view1/generate_eps
index 371e1d90ad31ae49101391d5a6f0e5c127e1193d..c863fbaddfb46cc5d1342d0ca9eb2766b31ac3dd 100755
--- a/scripts/view1/generate_eps
+++ b/scripts/view1/generate_eps
@@ -9,8 +9,8 @@ function print_the_help {
exit
}
-FILE_TAG="view1"
-INPUT_FILE="g4_0000.prim"
+FILE_TAG="view01"
+INPUT_FILE="../../g4_0000.prim"
POSITIONAL=()
@@ -44,6 +44,7 @@ done
set -- "${POSITIONAL[@]}" # restore positional parameters
+# Side view
dawncut 1 0 0 1 ${INPUT_FILE} ${FILE_TAG}_temp0.prim
dawncut -1 0 0 1 ${FILE_TAG}_temp0.prim ${FILE_TAG}.prim
dawn -d ${FILE_TAG}.prim
@@ -55,3 +56,15 @@ gs -o ${FILE_TAG}.pdf -sDEVICE=pdfwrite \
pdftoppm ${FILE_TAG}.pdf ${FILE_TAG} -png -singlefile -cropbox
+# Top view
+dawncut 0 1 0 1 ${INPUT_FILE} ${FILE_TAG}_temp0.prim
+dawncut 0 -1 0 1 ${FILE_TAG}_temp0.prim ${FILE_TAG}.prim
+../../bin/dawn_tweak --theta 270
+dawn -d ${FILE_TAG}.prim
+ps2pdf ${FILE_TAG}.eps ${FILE_TAG}_top_full.pdf
+gs -o ${FILE_TAG}_top.pdf -sDEVICE=pdfwrite \
+ -c "[/CropBox [51 250 550 590] /PAGES pdfmark" \
+ -f ${FILE_TAG}_top_full.pdf
+
+pdftoppm ${FILE_TAG}_top.pdf ${FILE_TAG}_top -png -singlefile -cropbox
+
diff --git a/src/ForwardRICH_geo.cpp b/src/ForwardRICH_geo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6520c30aac30492a8950c1590661bce57e089e47
--- /dev/null
+++ b/src/ForwardRICH_geo.cpp
@@ -0,0 +1,263 @@
+//==========================================================================
+// 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)
+