athena/src/DRICH_geo.cpp

0001 //==========================================================================
0002 //  dRICh: Dual Ring Imaging Cherenkov Detector
0003 //--------------------------------------------------------------------------
0004 //
0005 // Author: Christopher Dilks (Duke University)
0006 //
0007 // - Design Adapted from Standalone Fun4all and GEMC implementations
0008 //   [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
0009 //
0010 //==========================================================================
0011
0012 #include <XML/Helper.h>
0013 #include "TMath.h"
0014 #include "TString.h"
0015 #include "GeometryHelpers.h"
0016 #include "Math/Point2D.h"
0017 #include "DDRec/Surface.h"
0018 #include "DDRec/DetectorData.h"
0019 #include "DD4hep/OpticalSurfaces.h"
0020 #include "DD4hep/DetFactoryHelper.h"
0021 #include "DD4hep/Printout.h"
0022
0023 using namespace dd4hep;
0024 using namespace dd4hep::rec;
0025
0026 // create the detector
0027 static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens) {
0028   xml::DetElement detElem = handle;
0029   std::string detName = detElem.nameStr();
0030   int detID = detElem.id();
0031
0032   DetElement det(detName, detID);
0033   xml::Component dims = detElem.dimensions();
0034   OpticalSurfaceManager surfMgr = desc.surfaceManager();
0035
0036   // attributes -----------------------------------------------------------
0037   // - vessel
0038   double  vesselZmin       =  dims.attr<double>(_Unicode(zmin));
0039   double  vesselLength     =  dims.attr<double>(_Unicode(length));
0040   double  vesselRmin0      =  dims.attr<double>(_Unicode(rmin0));
0041   double  vesselRmin1      =  dims.attr<double>(_Unicode(rmin1));
0042   double  vesselRmax0      =  dims.attr<double>(_Unicode(rmax0));
0043   double  vesselRmax1      =  dims.attr<double>(_Unicode(rmax1));
0044   double  vesselRmax2      =  dims.attr<double>(_Unicode(rmax2));
0045   double  snoutLength      =  dims.attr<double>(_Unicode(snout_length));
0046   int     nSectors         =  dims.attr<int>(_Unicode(nsectors));
0047   double  wallThickness    =  dims.attr<double>(_Unicode(wall_thickness));
0048   double  windowThickness  =  dims.attr<double>(_Unicode(window_thickness));
0049   auto    vesselMat        =  desc.material(detElem.attr<std::string>(_Unicode(material)));
0050   auto    gasvolMat        =  desc.material(detElem.attr<std::string>(_Unicode(gas)));
0051   auto    vesselVis        =  desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
0052   auto    gasvolVis        =  desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
0053   // - radiator (applies to aerogel and filter)
0054   auto    radiatorElem        =  detElem.child(_Unicode(radiator));
0055   double  radiatorRmin        =  radiatorElem.attr<double>(_Unicode(rmin));
0056   double  radiatorRmax        =  radiatorElem.attr<double>(_Unicode(rmax));
0057   double  radiatorPhiw        =  radiatorElem.attr<double>(_Unicode(phiw));
0058   double  radiatorPitch       =  radiatorElem.attr<double>(_Unicode(pitch));
0059   double  radiatorFrontplane  =  radiatorElem.attr<double>(_Unicode(frontplane));
0060   // - aerogel
0061   auto    aerogelElem       =  radiatorElem.child(_Unicode(aerogel));
0062   auto    aerogelMat        =  desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
0063   auto    aerogelVis        =  desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
0064   double  aerogelThickness  =  aerogelElem.attr<double>(_Unicode(thickness));
0065   // - filter
0066   auto    filterElem       =  radiatorElem.child(_Unicode(filter));
0067   auto    filterMat        =  desc.material(filterElem.attr<std::string>(_Unicode(material)));
0068   auto    filterVis        =  desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
0069   double  filterThickness  =  filterElem.attr<double>(_Unicode(thickness));
0070   // - mirror
0071   auto    mirrorElem       =  detElem.child(_Unicode(mirror));
0072   auto    mirrorMat        =  desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
0073   auto    mirrorVis        =  desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
0074   auto    mirrorSurf       =  surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
0075   double  mirrorBackplane  =  mirrorElem.attr<double>(_Unicode(backplane));
0076   double  mirrorThickness  =  mirrorElem.attr<double>(_Unicode(thickness));
0077   double  mirrorRmin       =  mirrorElem.attr<double>(_Unicode(rmin));
0078   double  mirrorRmax       =  mirrorElem.attr<double>(_Unicode(rmax));
0079   double  mirrorPhiw       =  mirrorElem.attr<double>(_Unicode(phiw));
0080   double  focusTuneZ       =  mirrorElem.attr<double>(_Unicode(focus_tune_z));
0081   double  focusTuneX       =  mirrorElem.attr<double>(_Unicode(focus_tune_x));
0082   // - sensor module
0083   auto    sensorElem       =  detElem.child(_Unicode(sensors)).child(_Unicode(module));
0084   auto    sensorMat        =  desc.material(sensorElem.attr<std::string>(_Unicode(material)));
0085   auto    sensorVis        =  desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
0086   auto    sensorSurf       =  surfMgr.opticalSurface(sensorElem.attr<std::string>(_Unicode(surface)));
0087   double  sensorSide       =  sensorElem.attr<double>(_Unicode(side));
0088   double  sensorGap        =  sensorElem.attr<double>(_Unicode(gap));
0089   double  sensorThickness  =  sensorElem.attr<double>(_Unicode(thickness));
0090   // - sensor sphere
0091   auto    sensorSphElem     =  detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
0092   double  sensorSphRadius   =  sensorSphElem.attr<double>(_Unicode(radius));
0093   double  sensorSphCenterX  =  sensorSphElem.attr<double>(_Unicode(centerx));
0094   double  sensorSphCenterZ  =  sensorSphElem.attr<double>(_Unicode(centerz));
0095   // - sensor sphere patch cuts
0096   auto    sensorSphPatchElem  =  detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
0097   double  sensorSphPatchPhiw  =  sensorSphPatchElem.attr<double>(_Unicode(phiw));
0098   double  sensorSphPatchRmin  =  sensorSphPatchElem.attr<double>(_Unicode(rmin));
0099   double  sensorSphPatchRmax  =  sensorSphPatchElem.attr<double>(_Unicode(rmax));
0100   double  sensorSphPatchZmin  =  sensorSphPatchElem.attr<double>(_Unicode(zmin));
0101   // - debugging switches
0102   int   debug_optics_mode  =  detElem.attr<int>(_Unicode(debug_optics));
0103   bool  debug_mirror       =  mirrorElem.attr<bool>(_Unicode(debug));
0104   bool  debug_sensors      =  sensorSphElem.attr<bool>(_Unicode(debug));
0105
0106   // if debugging optics, override some settings
0107   bool debug_optics = debug_optics_mode > 0;
0108   if(debug_optics) {
0109     printout(WARNING,"DRich_geo","DEBUGGING DRICH OPTICS");
0110     switch(debug_optics_mode) {
0111       case 1: vesselMat = aerogelMat = filterMat = sensorMat = gasvolMat = desc.material("VacuumOptical"); break;
0112       case 2: vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical"); break;
0113       default: printout(FATAL,"DRich_geo","UNKNOWN debug_optics_mode"); return det;
0114     };
0115     aerogelVis = sensorVis = mirrorVis;
0116     gasvolVis = vesselVis = desc.invisible();
0117   };
0118
0119
0120   // BUILD VESSEL ====================================================================
0121   /* - `vessel`: aluminum enclosure, the mother volume of the dRICh
0122    * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
0123    *   are children of `gasvol`
0124    * - the dRICh vessel geometry has two regions: the snout refers to the conic region
0125    *   in the front, housing the aerogel, while the tank refers to the cylindrical
0126    *   region, housing the rest of the detector components
0127    */
0128
0129   // derived attributes
0130   double tankLength = vesselLength - snoutLength;
0131   double vesselZmax = vesselZmin + vesselLength;
0132
0133   // snout solids
0134   double boreDelta = vesselRmin1 - vesselRmin0;
0135   double snoutDelta = vesselRmax1 - vesselRmax0;
0136   Cone vesselSnout(
0137       snoutLength/2.0,
0138       vesselRmin0,
0139       vesselRmax0,
0140       vesselRmin0 + boreDelta * snoutLength / vesselLength,
0141       vesselRmax1
0142       );
0143   Cone gasvolSnout(
0144       /* note: `gasvolSnout` extends a bit into the tank, so it touches `gasvolTank`
0145        * - the extension distance is equal to the tank `windowThickness`, so the
0146        *   length of `gasvolSnout` == length of `vesselSnout`
0147        * - the extension backplane radius is calculated using similar triangles
0148        */
0149       snoutLength/2.0,
0150       vesselRmin0 + wallThickness,
0151       vesselRmax0 - wallThickness,
0152       vesselRmin0 + boreDelta * (snoutLength-windowThickness) / vesselLength + wallThickness,
0153       vesselRmax1 - wallThickness + windowThickness * (vesselRmax1 - vesselRmax0) / snoutLength
0154       );
0155
0156   // tank solids
0157   Cone vesselTank(
0158       tankLength/2.0,
0159       vesselSnout.rMin2(),
0160       vesselRmax2,
0161       vesselRmin1,
0162       vesselRmax2
0163       );
0164   Cone gasvolTank(
0165       tankLength/2.0 - windowThickness,
0166       gasvolSnout.rMin2(),
0167       vesselRmax2 - wallThickness,
0168       vesselRmin1 + wallThickness,
0169       vesselRmax2 - wallThickness
0170       );
0171
0172   // snout + tank solids
0173   UnionSolid vesselUnion(
0174       vesselTank,
0175       vesselSnout,
0176       Position(0., 0., -vesselLength/2.)
0177       );
0178   UnionSolid gasvolUnion(
0179       gasvolTank,
0180       gasvolSnout,
0181       Position(0., 0., -vesselLength/2. + windowThickness)
0182       );
0183
0184   //  extra solids for `debug_optics` only
0185   Box vesselBox(1001,1001,1001);
0186   Box gasvolBox(1000,1000,1000);
0187
0188   // choose vessel and gasvol solids (depending on `debug_optics_mode` (0=disabled))
0189   Solid vesselSolid, gasvolSolid;
0190   switch(debug_optics_mode) {
0191     case 0:  vesselSolid=vesselUnion;  gasvolSolid=gasvolUnion;  break; // `!debug_optics`
0192     case 1:  vesselSolid=vesselBox;    gasvolSolid=gasvolBox;    break;
0193     case 2:  vesselSolid=vesselBox;    gasvolSolid=gasvolUnion;  break;
0194   };
0195
0196   // volumes
0197   Volume vesselVol(detName, vesselSolid, vesselMat);
0198   Volume gasvolVol(detName+"_gas", gasvolSolid, gasvolMat);
0199   vesselVol.setVisAttributes(vesselVis);
0200   gasvolVol.setVisAttributes(gasvolVis);
0201
0202   // reference positions
0203   // - the vessel is created such that the center of the cylindrical tank volume
0204   //   coincides with the origin; this is called the "origin position" of the vessel
0205   // - when the vessel (and its children volumes) is placed, it is translated in
0206   //   the z-direction to be in the proper ATHENA-integration location
0207   // - these reference positions are for the frontplane and backplane of the vessel,
0208   //   with respect to the vessel origin position
0209   auto originFront = Position(0., 0., -tankLength/2.0 - snoutLength );
0210   auto originBack =  Position(0., 0., tankLength/2.0 );
0211
0212   // initialize sensor centroids (used for mirror parameterization below); this is
0213   // the average (x,y,z) of the placed sensors, w.r.t. originFront
0214   double sensorCentroidX = 0;
0215   double sensorCentroidZ = 0;
0216   int sensorCount = 0;
0217
0218
0219   // sensitive detector type
0220   sens.setType("tracker");
0221
0222
0223   // BUILD RADIATOR ====================================================================
0224
0225   // attributes
0226   double airGap = 0.01*mm; // air gap between aerogel and filter (FIXME? actually it's currently a gas gap)
0227
0228   // solid and volume: create aerogel and filter
0229   Cone aerogelSolid(
0230       aerogelThickness/2,
0231       radiatorRmin,
0232       radiatorRmax,
0233       radiatorRmin + boreDelta  * aerogelThickness / vesselLength,
0234       radiatorRmax + snoutDelta * aerogelThickness / snoutLength
0235       );
0236   Cone filterSolid(
0237       filterThickness/2,
0238       radiatorRmin + boreDelta  * (aerogelThickness + airGap) / vesselLength,
0239       radiatorRmax + snoutDelta * (aerogelThickness + airGap) / snoutLength,
0240       radiatorRmin + boreDelta  * (aerogelThickness + airGap + filterThickness) / vesselLength,
0241       radiatorRmax + snoutDelta * (aerogelThickness + airGap + filterThickness) / snoutLength
0242       );
0243
0244   Volume aerogelVol( detName+"_aerogel", aerogelSolid, aerogelMat );
0245   Volume filterVol(  detName+"_filter",  filterSolid,  filterMat );
0246   aerogelVol.setVisAttributes(aerogelVis);
0247   filterVol.setVisAttributes(filterVis);
0248
0249   // aerogel placement and surface properties
0250   // TODO [low-priority]: define skin properties for aerogel and filter
0251   auto radiatorPos = Position(0., 0., radiatorFrontplane) + originFront;
0252   auto aerogelPV = gasvolVol.placeVolume(aerogelVol,
0253         Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
0254       * RotationY(radiatorPitch) // change polar angle to specified pitch // FIXME: probably broken (not yet in use anyway)
0255       );
0256   DetElement aerogelDE(det, "aerogel_de", 0);
0257   aerogelDE.setPlacement(aerogelPV);
0258   //SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
0259   //aerogelSkin.isValid();
0260
0261   // filter placement and surface properties
0262   if(!debug_optics) {
0263     auto filterPV = gasvolVol.placeVolume(filterVol,
0264           Translation3D(0., 0., airGap) // add an air gap
0265         * Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
0266         * RotationY(radiatorPitch) // change polar angle
0267         * Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
0268         );
0269     DetElement filterDE(det, "filter_de", 0);
0270     filterDE.setPlacement(filterPV);
0271     //SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
0272     //filterSkin.isValid();
0273   };
0274
0275
0276   // SECTOR LOOP //////////////////////////////////
0277   for(int isec=0; isec<nSectors; isec++) {
0278
0279     // debugging filters, limiting the number of sectors
0280     if( (debug_mirror||debug_sensors||debug_optics) && isec!=0) continue;
0281
0282     // sector rotation about z axis
0283     double sectorRotation = isec * 360/nSectors * degree;
0284     std::string secName = "sec" + std::to_string(isec);
0285
0286
0287
0288     // BUILD SENSORS ====================================================================
0289
0290     // if debugging sphere properties, restrict number of sensors drawn
0291     if(debug_sensors) { sensorSide = 2*M_PI*sensorSphRadius / 64; };
0292
0293     // solid and volume: single sensor module
0294     Box sensorSolid(sensorSide/2., sensorSide/2., sensorThickness/2.);
0295     Volume sensorVol(detName+"_sensor_"+secName, sensorSolid, sensorMat);
0296     sensorVol.setVisAttributes(sensorVis);
0297
0298     auto sensorSphPos = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;
0299
0300     // sensitivity
0301     if(!debug_optics) sensorVol.setSensitiveDetector(sens);
0302
0303     // SENSOR MODULE LOOP ------------------------
0304     /* ALGORITHM: generate sphere of positions
0305      * - NOTE: there are two coordinate systems here:
0306      *   - "global" the main ATHENA coordinate system
0307      *   - "generator" (vars end in `Gen`) is a local coordinate system for
0308      *     generating points on a sphere; it is related to the global system by
0309      *     a rotation; we do this so the "patch" (subset of generated
0310      *     positions) of sensors we choose to build is near the equator, where
0311      *     point distribution is more uniform
0312      * - PROCEDURE: loop over `thetaGen`, with subloop over `phiGen`, each divided evenly
0313      *   - the number of points to generate depends how many sensors (+`sensorGap`)
0314      *     can fit within each ring of constant `thetaGen` or `phiGen`
0315      *   - we divide the relevant circumference by the sensor
0316      *     size(+`sensorGap`), and this number is allowed to be a fraction,
0317      *     because likely we don't care about generating a full sphere and
0318      *     don't mind a "seam" at the overlap point
0319      *   - if we pick a patch of the sphere near the equator, and not near
0320      *     the poles or seam, the sensor distribution will appear uniform
0321      */
0322
0323     // initialize module number for this sector
0324     int imod=0;
0325
0326     // thetaGen loop: iterate less than "0.5 circumference / sensor size" times
0327     double nTheta = M_PI*sensorSphRadius / (sensorSide+sensorGap);
0328     for(int t=0; t<(int)(nTheta+0.5); t++) {
0329       double thetaGen = t/((double)nTheta) * M_PI;
0330
0331       // phiGen loop: iterate less than "circumference at this latitude / sensor size" times
0332       double nPhi = 2*M_PI * sensorSphRadius * std::sin(thetaGen) / (sensorSide+sensorGap);
0333       for(int p=0; p<(int)(nPhi+0.5); p++) {
0334         double phiGen = p/((double)nPhi) * 2*M_PI - M_PI; // shift to [-pi,pi]
0335
0336         // determine global phi and theta
0337         // - convert {radius,thetaGen,phiGen} -> {xGen,yGen,zGen}
0338         double xGen = sensorSphRadius * std::sin(thetaGen) * std::cos(phiGen);
0339         double yGen = sensorSphRadius * std::sin(thetaGen) * std::sin(phiGen);
0340         double zGen = sensorSphRadius * std::cos(thetaGen);
0341         // - convert {xGen,yGen,zGen} -> global {x,y,z} via rotation
0342         double x = zGen;
0343         double y = xGen;
0344         double z = yGen;
0345         // - convert global {x,y,z} -> global {phi,theta}
0346         double phi = std::atan2(y,x);
0347         double theta = std::acos(z/sensorSphRadius);
0348
0349         // shift global coordinates so we can apply spherical patch cuts
0350         double zCheck = z + sensorSphCenterZ;
0351         double xCheck = x + sensorSphCenterX;
0352         double yCheck = y;
0353         double rCheck = std::hypot(xCheck,yCheck);
0354         double phiCheck = std::atan2(yCheck,xCheck);
0355
0356         // patch cut
0357         bool patchCut =
0358           std::fabs(phiCheck) < sensorSphPatchPhiw
0359           && zCheck > sensorSphPatchZmin
0360           && rCheck > sensorSphPatchRmin
0361           && rCheck < sensorSphPatchRmax;
0362         if(debug_sensors) patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
0363         if(patchCut) {
0364
0365           // append sensor position to centroid calculation
0366           if(isec==0) {
0367             sensorCentroidX += xCheck;
0368             sensorCentroidZ += zCheck;
0369             sensorCount++;
0370           };
0371
0372           // placement (note: transformations are in reverse order)
0373           // - transformations operate on global coordinates; the corresponding
0374           //   generator coordinates are provided in the comments
0375           auto sensorPV = gasvolVol.placeVolume(sensorVol,
0376                 RotationZ(sectorRotation) // rotate about beam axis to sector
0377               * Translation3D(sensorSphPos.x(), sensorSphPos.y(), sensorSphPos.z()) // move sphere to reference position
0378               * RotationX(phiGen) // rotate about `zGen`
0379               * RotationZ(thetaGen) // rotate about `yGen`
0380               * Translation3D(sensorSphRadius, 0., 0.) // push radially to spherical surface
0381               * RotationY(M_PI/2) // rotate sensor to be compatible with generator coords
0382               * RotationZ(-M_PI/2) // correction for readout segmentation mapping
0383               );
0384
0385           // generate LUT for module number -> sensor position, for readout mapping tests
0386           //if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
0387
0388           // properties
0389           sensorPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);
0390           DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), (imod<<3)|isec ); // id must match IRTAlgorithm usage
0391           sensorDE.setPlacement(sensorPV);
0392           if(!debug_optics) {
0393             SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVol);
0394             sensorSkin.isValid();
0395           };
0396
0397           // increment sensor module number
0398           imod++;
0399
0400         }; // end patch cuts
0401       }; // end phiGen loop
0402     }; // end thetaGen loop
0403
0404     // calculate centroid sensor position
0405     if(isec==0) {
0406       sensorCentroidX /= sensorCount;
0407       sensorCentroidZ /= sensorCount;
0408     };
0409
0410     // END SENSOR MODULE LOOP ------------------------
0411
0412
0413     // BUILD MIRRORS ====================================================================
0414
0415     // derive spherical mirror parameters `(zM,xM,rM)`, for given image point
0416     // coordinates `(zI,xI)` and `dO`, defined as the z-distance between the
0417     // object and the mirror surface
0418     // - all coordinates are specified w.r.t. the object point coordinates
0419     // - this is point-to-point focusing, but it can be used to effectively steer
0420     //   parallel-to-point focusing
0421     double zM,xM,rM;
0422     auto FocusMirror = [&zM,&xM,&rM](double zI, double xI, double dO) {
0423       zM = dO*zI / (2*dO-zI);
0424       xM = dO*xI / (2*dO-zI);
0425       rM = dO - zM;
0426     };
0427
0428     // attributes, re-defined w.r.t. IP, needed for mirror positioning
0429     double zS = sensorSphCenterZ + vesselZmin; // sensor sphere attributes
0430     double xS = sensorSphCenterX;
0431     double rS = sensorSphRadius;
0432     double B = vesselZmax - mirrorBackplane; // distance between IP and mirror back plane
0433
0434     // focus 1: set mirror to focus IP on center of sensor sphere `(zS,xS)`
0435     /*double zF = zS;
0436     double xF = xS;
0437     FocusMirror(zF,xF,B);*/
0438
0439     // focus 2: move focal region along sensor sphere radius, according to `focusTuneLong`
0440     // - specifically, along the radial line which passes through the approximate centroid
0441     //   of the sensor region `(sensorCentroidZ,sensorCentroidX)`
0442     // - `focusTuneLong` is the distance to move, given as a fraction of `sensorSphRadius`
0443     // - `focusTuneLong==0` means `(zF,xF)==(zS,xS)`
0444     // - `focusTuneLong==1` means `(zF,xF)` will be on the sensor sphere, near the centroid
0445     /*
0446     double zC = sensorCentroidZ + vesselZmin;
0447     double xC = sensorCentroidX;
0448     double slopeF = (xC-xS) / (zC-zS);
0449     double thetaF = std::atan(std::fabs(slopeF));
0450     double zF = zS + focusTuneLong * sensorSphRadius * std::cos(thetaF);
0451     double xF = xS - focusTuneLong * sensorSphRadius * std::sin(thetaF);
0452     //FocusMirror(zF,xF,B);
0453
0454     // focus 3: move along line perpendicular to focus 2's radial line,
0455     // according to `focusTunePerp`, with the same numerical scale as `focusTuneLong`
0456     zF += focusTunePerp * sensorSphRadius * std::cos(M_PI/2-thetaF);
0457     xF += focusTunePerp * sensorSphRadius * std::sin(M_PI/2-thetaF);
0458     FocusMirror(zF,xF,B);
0459     */
0460
0461     // focus 4: use (z,x) coordinates for tune parameters
0462     double zF = zS + focusTuneZ;
0463     double xF = xS + focusTuneX;
0464     FocusMirror(zF,xF,B);
0465
0466     // re-define mirror attributes to be w.r.t vessel front plane
0467     double mirrorCenterZ = zM - vesselZmin;
0468     double mirrorCenterX = xM;
0469     double mirrorRadius = rM;
0470
0471     // spherical mirror patch cuts and rotation
0472     double mirrorThetaRot = std::asin(mirrorCenterX/mirrorRadius);
0473     double mirrorTheta1 = mirrorThetaRot - std::asin((mirrorCenterX-mirrorRmin)/mirrorRadius);
0474     double mirrorTheta2 = mirrorThetaRot + std::asin((mirrorRmax-mirrorCenterX)/mirrorRadius);
0475
0476     // if debugging, draw full sphere
0477     if(debug_mirror) { mirrorTheta1=0; mirrorTheta2=M_PI; /*mirrorPhiw=2*M_PI;*/ };
0478
0479     // solid : create sphere at origin, with specified angular limits;
0480     // phi limits are increased to fill gaps (overlaps are cut away later)
0481     Sphere mirrorSolid1(
0482         mirrorRadius,
0483         mirrorRadius + mirrorThickness,
0484         mirrorTheta1,
0485         mirrorTheta2,
0486         -40*degree,
0487         40*degree
0488         );
0489
0490     /* CAUTION: if any of the relative placements or boolean operations below
0491      * are changed, you MUST make sure this does not break access to the sphere
0492      * primitive and positioning in Juggler `IRTAlgorithm`; cross check the
0493      * mirror sphere attributes carefully!
0494      */
0495     /*
0496     // PRINT MIRROR ATTRIBUTES (before any sector z-rotation)
0497     printf("zM = %f\n",zM); // sphere centerZ, w.r.t. IP
0498     printf("xM = %f\n",xM); // sphere centerX, w.r.t. IP
0499     printf("rM = %f\n",rM); // sphere radius
0500     */
0501
0502     // mirror placement transformation (note: transformations are in reverse order)
0503     auto mirrorPos = Position(mirrorCenterX, 0., mirrorCenterZ) + originFront;
0504     Transform3D mirrorPlacement(
0505           Translation3D(mirrorPos.x(), mirrorPos.y(), mirrorPos.z()) // re-center to specified position
0506         * RotationY(-mirrorThetaRot) // rotate about vertical axis, to be within vessel radial walls
0507         );
0508
0509     // cut overlaps with other sectors using "pie slice" wedges, to the extent specified
0510     // by `mirrorPhiw`
0511     Tube pieSlice( 0.01*cm, vesselRmax2, tankLength/2.0, -mirrorPhiw/2.0, mirrorPhiw/2.0);
0512     IntersectionSolid mirrorSolid2( pieSlice, mirrorSolid1, mirrorPlacement );
0513
0514     // mirror volume, attributes, and placement
0515     Volume mirrorVol(detName+"_mirror_"+secName, mirrorSolid2, mirrorMat);
0516     mirrorVol.setVisAttributes(mirrorVis);
0517     auto mirrorPV2 = gasvolVol.placeVolume(mirrorVol,
0518           RotationZ(sectorRotation) // rotate about beam axis to sector
0519         * Translation3D(0,0,0)
0520         );
0521
0522     // properties
0523     DetElement mirrorDE(det, Form("mirror_de%d", isec), isec);
0524     mirrorDE.setPlacement(mirrorPV2);
0525     SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
0526     mirrorSkin.isValid();
0527
0528
0529   }; // END SECTOR LOOP //////////////////////////
0530
0531
0532   // place gas volume
0533   PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol,Position(0, 0, 0));
0534   DetElement gasvolDE(det, "gasvol_de", 0);
0535   gasvolDE.setPlacement(gasvolPV);
0536
0537   // place mother volume (vessel)
0538   Volume motherVol = desc.pickMotherVolume(det);
0539   PlacedVolume vesselPV = motherVol.placeVolume(vesselVol,
0540       Position(0, 0, vesselZmin) - originFront
0541       );
0542   vesselPV.addPhysVolID("system", detID);
0543   det.setPlacement(vesselPV);
0544
0545   return det;
0546 };
0547
0548 // clang-format off
0549 DECLARE_DETELEMENT(athena_DRICH, createDetector)