athena/src/GaseousRICH_geo.cpp

0001 //==========================================================================
0002 //  Gaseous Ring Imaging Cherenkov Detector
0003 //--------------------------------------------------------------------------
0004 //
0005 // Author: C. Peng (ANL)
0006 // Date: 09/30/2020
0007 //
0008 //==========================================================================
0009
0010 #include <XML/Helper.h>
0011 #include "TMath.h"
0012 #include "TString.h"
0013 #include "GeometryHelpers.h"
0014 #include "Math/Point2D.h"
0015 #include "DDRec/Surface.h"
0016 #include "DDRec/DetectorData.h"
0017 #include "DD4hep/OpticalSurfaces.h"
0018 #include "DD4hep/DetFactoryHelper.h"
0019 #include "DD4hep/Printout.h"
0020
0021 using namespace std;
0022 using namespace dd4hep;
0023 using namespace dd4hep::rec;
0024
0025 // headers
0026 void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset);
0027 void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset);
0028 void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens);
0029
0030 // helper function to get x, y, z if defined in a xml component
0031 template<class XmlComp>
0032 Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
0033 {
0034     Position pos(0., 0., 0.);
0035     if (comp.hasChild(name)) {
0036         auto child = comp.child(name);
0037         pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
0038         pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
0039         pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
0040     }
0041     return pos;
0042 }
0043
0044 // create the detector
0045 static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
0046 {
0047     xml::DetElement detElem = handle;
0048     std::string detName = detElem.nameStr();
0049     int detID = detElem.id();
0050
0051     DetElement det(detName, detID);
0052     xml::Component dims = detElem.dimensions();
0053
0054     // build a big envelope for all the components, filled with optical material to ensure transport of optical photons
0055     double z0 = dims.z0();
0056     double length = dims.length();
0057     double rmin = dims.rmin();
0058     double rmax0 = dims.attr<double>(_Unicode(rmax0));
0059     double rmax1 = dims.attr<double>(_Unicode(rmax1));
0060     double rmax2 = dims.attr<double>(_Unicode(rmax2));
0061     double snout_length = dims.attr<double>(_Unicode(snout_length));
0062     // fill envelope with radiator materials (need optical property for optical photons)
0063     auto gasMat = desc.material(dd4hep::getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
0064
0065     Cone snout(snout_length/2.0, rmin, rmax0, rmin, rmax1);
0066     Tube tank(rmin, rmax2, (length - snout_length)/2., 0., 2*M_PI);
0067     // shift the snout to the left side of tank
0068     UnionSolid envShape(tank, snout, Position(0., 0., -length/2.));
0069     // some reference points
0070     auto snout_front = Position(0., 0., -(length + snout_length)/2.);
0071     // auto snout_end = Position(0., 0., -(length - snout_length)/2.);  // tank_front
0072     // auto tank_end = Position(0., 0., (length - snout_length)/2.);
0073
0074     Volume envVol(detName + "_envelope", envShape, gasMat);
0075     envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
0076
0077     // sensitive detector type
0078     sens.setType("tracker");
0079
0080     // @TODO: place a radiator
0081     // build_radiator(desc, envVol, detElement.child(_Unicode(radiator)), snout_front);
0082
0083     // place mirrors
0084     build_mirrors(desc, det, envVol, detElem.child(_Unicode(mirrors)), snout_front);
0085
0086     // place photo-sensors
0087     build_sensors(desc, envVol, detElem.child(_Unicode(sensors)), snout_front, sens);
0088
0089     Volume motherVol = desc.pickMotherVolume(det);
0090     PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0) - snout_front);
0091     envPV.addPhysVolID("system", detID);
0092     det.setPlacement(envPV);
0093
0094     return det;
0095 }
0096
0097 // @TODO: implement a radiator, now the envelope serves as the radiator
0098 void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset)
0099 {
0100     // place holder
0101 }
0102
0103 // place mirrors
0104 void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset)
0105 {
0106     double thickness = dd4hep::getAttrOrDefault<double>(plm, _Unicode(dz), 1.*dd4hep::mm);
0107     auto mat = desc.material(plm.attr<std::string>(_Unicode(material)));
0108     auto vis = desc.visAttributes(plm.attr<std::string>(_Unicode(vis)));
0109
0110     // optical surface
0111     OpticalSurfaceManager surfMgr = desc.surfaceManager();
0112     auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(plm, _Unicode(surface), "MirrorOpticalSurface"));
0113
0114     // placements
0115     auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
0116     auto grot = get_xml_xyz(plm, _Unicode(position));
0117     int imod = 1;
0118     for (xml::Collection_t mir(plm, _Unicode(mirror)); mir; ++mir, ++imod) {
0119         double rmin = mir.attr<double>(_Unicode(rmin));
0120         double rmax = mir.attr<double>(_Unicode(rmax));
0121         double phiw = dd4hep::getAttrOrDefault<double>(mir, _Unicode(phiw), 2.*M_PI);
0122         double rotz = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotz), 0.);
0123         double roty = dd4hep::getAttrOrDefault<double>(mir, _Unicode(roty), 0.);
0124         double rotx = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotx), 0.);
0125
0126         Volume vol(Form("mirror_v_dummy%d", imod));
0127         vol.setMaterial(mat);
0128         vol.setVisAttributes(vis);
0129         // mirror curvature
0130         double curve = dd4hep::getAttrOrDefault<double>(mir, _Unicode(curve), 0.);
0131         // spherical mirror
0132         if (curve > 0.) {
0133             double th1 = std::asin(rmin/curve);
0134             double th2 = std::asin(rmax/curve);
0135             vol.setSolid(Sphere(curve, curve + thickness, th1, th2, 0., phiw));
0136         // plane mirror
0137         } else {
0138             vol.setSolid(Tube(rmin, rmax, thickness/2., 0., phiw));
0139         }
0140         // transforms are in a reverse order
0141         Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
0142                        * RotationZYX(grot.z(), grot.y(), grot.x())
0143                        * RotationZ(rotz)                // rotation of the piece
0144                        * RotationY(roty)                // rotation of the piece
0145                        * RotationX(rotx)                // rotation of the piece
0146                        * Translation3D(0., 0., -curve)  // move spherical shell to origin (no move for planes)
0147                        * RotationZ(-phiw/2.);           // center phi angle to 0. (-phiw/2., phiw/2.)
0148         auto pv = env.placeVolume(vol, tr);
0149         DetElement de(sdet, Form("mirror_de%d", imod), imod);
0150         de.setPlacement(pv);
0151         SkinSurface skin(desc, de, Form("mirror_optical_surface%d", imod), surf, vol);
0152         skin.isValid();
0153     }
0154 }
0155
0156 // place photo-sensors
0157 void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens)
0158 {
0159     // build sensor unit geometry
0160     auto mod = plm.child(_Unicode(module));
0161     double sx = mod.attr<double>(_Unicode(sx));
0162     double sy = mod.attr<double>(_Unicode(sy));
0163     double sz = mod.attr<double>(_Unicode(sz));
0164     double gap = mod.attr<double>(_Unicode(gap));
0165     auto mat = desc.material(mod.attr<std::string>(_Unicode(material)));
0166     auto vis = desc.visAttributes(mod.attr<std::string>(_Unicode(vis)));
0167
0168     Box sensor(sx/2., sy/2., sz/2.);
0169     Volume svol("sensor_v", sensor, mat);
0170     svol.setVisAttributes(vis);
0171
0172     // a thin layer of cherenkov gas for accepting optical photons
0173     auto opt = plm.child(_Unicode(optical));
0174     double opthick = opt.attr<double>(_Unicode(thickness));
0175     auto opmat = desc.material(opt.attr<std::string>(_Unicode(material)));
0176
0177     Box opshape(sx/2., sy/2., sz/2. + opthick/2.);
0178     Volume opvol("sensor_v_optical", opshape, opmat);
0179     opvol.placeVolume(svol, Position(0., 0., 0.));
0180     opvol.setSensitiveDetector(sens);
0181
0182     // photo-detector plane envelope
0183     auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
0184     auto grot = get_xml_xyz(plm, _Unicode(position));
0185     int isec = 1;
0186     for (xml::Collection_t sec(plm, _Unicode(sector)); sec; ++sec, ++isec) {
0187         double rmin = sec.attr<double>(_Unicode(rmin));
0188         double rmax = sec.attr<double>(_Unicode(rmax));
0189         double phiw = dd4hep::getAttrOrDefault<double>(sec, _Unicode(phiw), 2.*M_PI);
0190         double rotz = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotz), 0.);
0191         double roty = dd4hep::getAttrOrDefault<double>(sec, _Unicode(roty), 0.);
0192         double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
0193
0194         // fill sensors to the piece
0195         auto points = athena::geo::fillRectangles({0., 0.}, sx + gap, sy + gap, rmin - gap, rmax + gap, -phiw/2., phiw/2.);
0196         int imod = 1;
0197         for (auto &p : points) {
0198             // transofrms are in a reversed order
0199             Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
0200                            * RotationZYX(grot.z(), grot.y(), grot.x())
0201                            * RotationZ(rotz)                    // rotation of the sector
0202                            * RotationY(roty)                    // rotation of the sector
0203                            * RotationX(rotx)                    // rotation of the sector
0204                            * Translation3D(p.x(), p.y(), 0.);   // place modules in each sector
0205             auto pv = env.placeVolume(opvol, tr);
0206             pv.addPhysVolID("sector", isec).addPhysVolID("module", imod++);
0207         }
0208     }
0209 }
0210
0211 //@}
0212
0213 // clang-format off
0214 DECLARE_DETELEMENT(athena_GaseousRICH, createDetector)
0215