ForwardRICH_geo.cpp 10.2 KB
Newer Older
Chao Peng's avatar
Chao Peng committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
//==========================================================================
//  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
Chao Peng's avatar
Chao Peng committed
28
inline bool in_ring(const Point &pt, double side, double rmin, double rmax, double phmin, double phmax)
Chao Peng's avatar
Chao Peng committed
29
{
Chao Peng's avatar
Chao Peng committed
30
31
32
33
    if (pt.r() > rmax || pt.r() < rmin) {
        return false;
    }

Chao Peng's avatar
Chao Peng committed
34
35
36
37
38
39
40
41
    // 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) {
Chao Peng's avatar
Chao Peng committed
42
        if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
Chao Peng's avatar
Chao Peng committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
            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
Chao Peng's avatar
Chao Peng committed
62
63
void add_square(Point p, std::vector<Point> &res, double lside, double rmin, double rmax,
                double phmin, double phmax)
Chao Peng's avatar
Chao Peng committed
64
65
{
    // outside of the ring or overlapping
Chao Peng's avatar
Chao Peng committed
66
    if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, res)) {
Chao Peng's avatar
Chao Peng committed
67
68
69
70
71
72
        return;
    }

    res.emplace_back(p);

    // check adjacent squares
Chao Peng's avatar
Chao Peng committed
73
74
75
76
    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);
Chao Peng's avatar
Chao Peng committed
77
78
79
}

// fill squares
Chao Peng's avatar
Chao Peng committed
80
81
std::vector<Point> fill_squares(Point ref, double lside, double rmin, double rmax,
                                double phmin = 0., double phmax = 2.*M_PI)
Chao Peng's avatar
Chao Peng committed
82
83
84
85
86
{
    // 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);

Chao Peng's avatar
Chao Peng committed
87
    auto find_seed = [] (const Point &ref, int n, double side, double rmin, double rmax, double phmin, double phmax) {
Chao Peng's avatar
Chao Peng committed
88
89
90
        for (int ix = -n; ix < n; ++ix) {
            for (int iy = -n; iy < n; ++iy) {
                Point pt(ref.x() + ix*side, ref.y() + iy*side);
Chao Peng's avatar
Chao Peng committed
91
                if (in_ring(pt, side, rmin, rmax, phmin, phmax)) {
Chao Peng's avatar
Chao Peng committed
92
93
94
95
96
97
98
99
                    return pt;
                }
            }
        }
        return ref;
    };

    std::vector<Point> res;
Chao Peng's avatar
Chao Peng committed
100
101
    ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax, phmin, phmax);
    add_square(ref, res, lside, rmin, rmax, phmin, phmax);
Chao Peng's avatar
Chao Peng committed
102
103
104
105
106
107
108
109
110
111
112
113
114
    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();
Chao Peng's avatar
Chao Peng committed
115
    xml::Component rads = detElem.child(_Unicode(radiator));
Chao Peng's avatar
Chao Peng committed
116
117
118
119
120
    xml::Component mir = detElem.child(_Unicode(mirror));
    xml::Component mcp = detElem.child(_Unicode(mcppmt));

    // dimensions
    double z0 = dims.z0();
Chao Peng's avatar
Chao Peng committed
121
122
123
124
    double length = dims.length();
    double rmin = dims.rmin();
    double rmax1 = dims.attr<double>(_Unicode(rmax1));
    double rmax2 = dims.attr<double>(_Unicode(rmax2));
Chao Peng's avatar
Chao Peng committed
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

    // 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());
Chao Peng's avatar
Chao Peng committed
141
    auto gasMat = desc.material(rads.materialStr());
Chao Peng's avatar
Chao Peng committed
142
143
    auto mcpMat = desc.material(mcp.materialStr());

Chao Peng's avatar
Chao Peng committed
144
145
    // constants
    auto richCenterAngle = std::atan((rmin + (rmax2 - rmin)/2.)/mirZ);
Whitney Armstrong's avatar
Whitney Armstrong committed
146
    //std::cout << richCenterAngle*180./M_PI << std::endl;
Chao Peng's avatar
Chao Peng committed
147

Chao Peng's avatar
Chao Peng committed
148
    // an envelope for the detector
Chao Peng's avatar
Chao Peng committed
149
150
    // use a complicated shape to avoid conflict with the other parts
    // cone for radiator and the first set of mirrors
Chao Peng's avatar
Chao Peng committed
151
152
153
154
155
    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);
Chao Peng's avatar
Chao Peng committed
156

Chao Peng's avatar
Chao Peng committed
157
    UnionSolid envShape(env1, env2, Position(0., 0., pZ));
Chao Peng's avatar
Chao Peng committed
158

Chao Peng's avatar
Chao Peng committed
159
160
    Volume envVol(detName + "_envelope", envShape, gasMat);
    envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
Chao Peng's avatar
Chao Peng committed
161
162

    // ---------------
Chao Peng's avatar
Chao Peng committed
163
    // spherical mirrors inside it
Chao Peng's avatar
Chao Peng committed
164
165
166
167
168
169
170
171
172
173
174
175
176
    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));
Chao Peng's avatar
Chao Peng committed
177
178
179
180
181
182
183
184
        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);
Chao Peng's avatar
Chao Peng committed
185
        mirVol.setVisAttributes(desc.visAttributes(mir.visStr()));
Chao Peng's avatar
Chao Peng committed
186
187
188
189
190
191
192
193
194
195
196
        // 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
Chao Peng's avatar
Chao Peng committed
197
                           * RotationX(180*degree)
Chao Peng's avatar
Chao Peng committed
198
199
200
201
202
203
204
205
206
                           * 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.)
Chao Peng's avatar
Chao Peng committed
207
            mirPV = envVol.placeVolume(mirVol, tr);
Chao Peng's avatar
Chao Peng committed
208
        }
Chao Peng's avatar
Chao Peng committed
209
        mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", imod);
Chao Peng's avatar
Chao Peng committed
210
        DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
Chao Peng's avatar
Chao Peng committed
211
212
213
214
215
216
217
        mirDE.setPlacement(mirPV);
        SkinSurface mirSurfBorder(desc, mirDE, Form("RICHmirror%d", imod), mirSurf, mirVol);
        mirSurfBorder.isValid();
    }
    ilayer++;

    // ---------------
Chao Peng's avatar
Chao Peng committed
218
    // photo-detector unit
Chao Peng's avatar
Chao Peng committed
219
220
221
222
223
224
225
226
227
228
229
230
231
    // 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);

Chao Peng's avatar
Chao Peng committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
    // 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);
Chao Peng's avatar
Chao Peng committed
251
252
253
254
255
256
257
258
259
260
261
262
263
    }
    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(ForwardRICH, createDetector)