MRich_geo.cpp

//
// Author     : Whit Armstrong (warmstrong@anl.gov)
//
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
#include "GeometryHelpers.h"
#include "Math/Vector3D.h"
#include "Math/AxisAngle.h"
#include "Math/VectorUtil.h"

using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;

using Placements = vector<PlacedVolume>;



static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDetector sens){
  xml_det_t      x_det    = e;
  Material       air      = description.material("AirOptical");
  Material       vacuum   = description.vacuum();
  string         det_name = x_det.nameStr();
  DetElement     sdet(det_name, x_det.id());
  Assembly       assembly(det_name);
  sens.setType("tracker");
  OpticalSurfaceManager surfMgr = description.surfaceManager();

  bool projective = getAttrOrDefault(x_det, _Unicode(projective), false);
  bool reflect    = x_det.reflect(true);

  PlacedVolume pv;

  map<string, Volume>     modules;
  map<string, Placements> sensitives;
  map<string, Volume>     module_assemblies;
  std::map<std::string,DetElement> module_assembly_delements;

  int                     n_sensor = 1;

  // dimensions
  xml::Component dims   = x_det.dimensions();
  auto           rmin   = dims.rmin();
  auto           rmax   = dims.rmax();
  auto           length = dims.length();
  auto           zmin   = dims.zmin();
  auto           zpos   = zmin + length / 2;

  // envelope
  Tube   envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
  Volume envVol("MRICH_Envelope", envShape, air);
  envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
  if (x_det.hasChild(_Unicode(envelope))) {
    xml_comp_t x_envelope = x_det.child(_Unicode(envelope));
    double thickness  = x_envelope.thickness();
    Material material = description.material(x_envelope.materialStr());
    Tube   envInsideShape(rmin + thickness, rmax - thickness, length / 2. - thickness);
    SubtractionSolid envShellShape(envShape, envInsideShape);
    Volume envShell("MRICH_Envelope_Inside", envShellShape, material);
    envVol.placeVolume(envShell);
  }

  // expect only one module (for now)
  xml_comp_t x_mod = x_det.child(_U(module));
  string     mod_name            = x_mod.nameStr();
  double     mod_width           = getAttrOrDefault(x_mod, _U(width), 130.0 * mm);
  double     mod_height          = getAttrOrDefault(x_mod, _U(height), 130.0 * mm);
  double     mod_length          = getAttrOrDefault(x_mod, _U(length), 130.0 * mm);

  // module
  Box    m_solid(mod_width / 2.0, mod_height / 2.0, mod_length / 2.0);
  Volume m_volume(mod_name, m_solid, air);
  m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
  DetElement mod_de( mod_name + std::string("_mod_") + std::to_string(1), 1);
  double z_placement = - mod_length / 2.0;

  // todo module frame
  if (x_mod.hasChild(_Unicode(frame))) {
    xml_comp_t x_frame    = x_mod.child(_Unicode(frame));
    double     frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
    Box        frame_inside(mod_width / 2.0 - frame_thickness, mod_height / 2.0 - frame_thickness, mod_length / 2.0 - frame_thickness);
    SubtractionSolid frame_solid(m_solid, frame_inside);
    Material   frame_mat       = description.material(x_frame.materialStr());
    Volume     frame_vol(mod_name+"_frame", frame_solid, frame_mat);
    auto       frame_vis       = getAttrOrDefault<std::string>(x_frame, _U(vis), std::string("GrayVis"));
    frame_vol.setVisAttributes(description.visAttributes(frame_vis));
    // update position
    z_placement += frame_thickness / 2.0;
    // place volume
    m_volume.placeVolume(frame_vol);
    // update position
    z_placement += frame_thickness / 2.0;
  }

  // aerogel box
  if (x_mod.hasChild(_Unicode(aerogel))) {
    xml_comp_t x_aerogel  = x_mod.child(_Unicode(aerogel));
    double     aerogel_width   = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
    double     aerogel_length  = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
    Material   aerogel_mat     = description.material(x_aerogel.materialStr());
    auto       aerogel_vis     = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));

    xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
    double     foam_thickness  = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
    Material   foam_mat        = description.material(x_aerogel_frame.materialStr());
    auto       foam_vis        = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));

    // foam frame
    Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0);
    Box foam_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0);
    SubtractionSolid foam_frame_solid(foam_box, foam_sub_box, Position(0, 0, foam_thickness));
    Volume           foam_vol(mod_name+"_aerogel_frame", foam_frame_solid, foam_mat);
    foam_vol.setVisAttributes(description.visAttributes(foam_vis));

    // aerogel
    Box              aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0);
    Volume           aerogel_vol(mod_name+"_aerogel", aerogel_box, aerogel_mat);
    aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis));

    // update position
    z_placement += (aerogel_length + foam_thickness) / 2.0;
    // place foam frame
    pv = m_volume.placeVolume(foam_vol,Position(0,0,z_placement));
    // place aerogel
    z_placement += foam_thickness / 2.0;
    pv = m_volume.placeVolume(aerogel_vol,Position(0,0,z_placement));
    DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
    aerogel_de.setPlacement(pv);
    // update position
    z_placement += aerogel_length / 2.0;

    // optical surfaces
    auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
    SkinSurface skin_surf(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol);
    skin_surf.isValid();
  }

  // Fresnel Lens
  if (x_mod.hasChild(_Unicode(lens))) {
    xml_comp_t x_lens     = x_mod.child(_Unicode(lens));

    //  - The lens has a constant groove pitch (delta r) as opposed to fixing the groove height.
    //  - The lens area outside of the effective diamtere is flat.
    //  - The grooves are not curved, rather they are polycone shaped, ie a flat approximating the curvature.
    auto   lens_vis       = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue"));
    double groove_pitch   = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm);// 0.5 * mm);
    double lens_f         = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0*2.54*cm);
    double eff_diameter   = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * mm);
    double lens_width     = getAttrOrDefault(x_lens, _Unicode(width), 6.7*2.54*cm);
    double center_thickness = getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm);//2.0 * mm);

    double n_acrylic        = 1.49;
    double lens_curvature   = 1.0 / (lens_f*(n_acrylic - 1.0)); //confirmed
    double full_ring_rmax   = std::min(eff_diameter / 2.0, lens_width/2.0);

    double N_grooves        = std::ceil((full_ring_rmax) / groove_pitch);
    double groove_last_rmin = (N_grooves - 1) * groove_pitch;
    double groove_last_rmax = N_grooves * groove_pitch;

    auto   groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); };
    double lens_thickness = groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness;

    Material         lens_mat = description.material(x_lens.materialStr());
    Box              lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0);
    SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness));

    Assembly lens_vol(mod_name + "_lens");
    Volume   flatpart_lens_vol( "flatpart_lens", flat_lens, lens_mat);
    lens_vol.placeVolume(flatpart_lens_vol);

    int    i_groove           = 0;
    double groove_rmax        = groove_pitch;
    double groove_rmin        = 0;

    while ( groove_rmax <= full_ring_rmax ) {
      double   dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin);
      Polycone groove_solid(0, 2.0 * M_PI,
                            {groove_rmin, groove_rmin, groove_rmin},
                            {groove_rmax, groove_rmax, groove_rmin},
                            {-lens_thickness/2.0, lens_thickness/2.0-dZ, lens_thickness/2.0});
      Volume   lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat);
      lens_vol.placeVolume(lens_groove_vol);

      i_groove++;
      groove_rmin = (i_groove  )*groove_pitch;
      groove_rmax = (i_groove+1)*groove_pitch;
    }

    lens_vol.setVisAttributes(description.visAttributes(lens_vis));

    // update position
    z_placement += lens_thickness/2.0;
    // place volume
    pv = m_volume.placeVolume(lens_vol,Position(0,0,z_placement));
    DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
    lens_de.setPlacement(pv);
    // update position
    z_placement += lens_thickness/2.0;

    // optical surfaces
    auto lens_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
    SkinSurface skin_surf(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), lens_surf, lens_vol);
    skin_surf.isValid();
  }

  // mirror
  if (x_mod.hasChild(_Unicode(space))) {
    xml_comp_t x_space = x_mod.child(_Unicode(space));
    z_placement += getAttrOrDefault(x_space, _U(thickness), 0.0 * mm);
  }

  // mirror
  if (x_mod.hasChild(_Unicode(mirror))) {
    xml_comp_t x_mirror   = x_mod.child(_Unicode(mirror));
    auto   mirror_vis     = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray"));
    double mirror_x1      = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm);
    double mirror_x2      = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm);
    double mirror_length  = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm);
    double mirror_thickness  = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm);
    double outer_x1 = (mirror_x1+mirror_thickness)/2.0;
    double outer_x2 = (mirror_x2+mirror_thickness)/2.0;
    Trd2   outer_mirror_trd(outer_x1, outer_x2, outer_x1,  outer_x2, mirror_length/2.0);
    Trd2   inner_mirror_trd(mirror_x1 / 2.0,  mirror_x2 / 2.0, mirror_x1 / 2.0,mirror_x2 / 2.0, mirror_length/2.0+0.1*mm);
    SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd);
    Material mirror_mat        = description.material(x_mirror.materialStr());
    Volume   mirror_vol(mod_name+"_mirror", mirror_solid, mirror_mat);

    // update position
    z_placement += mirror_length/2.0;
    // place volume
    pv = m_volume.placeVolume(mirror_vol,Position(0,0,z_placement));
    DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
    mirror_de.setPlacement(pv);
    // update position
    z_placement += mirror_length/2.0;

    // optical surfaces
    auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
    SkinSurface skin_surf(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol);
    skin_surf.isValid();
  }

  // photon detector
  if (x_mod.hasChild(_Unicode(photodet))) {
    xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
    auto       photodet_vis       = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed"));
    double     photodet_width     = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm);
    double     photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm);
    Material   photodet_mat       = description.material(x_photodet.materialStr());
    Box        window_box(photodet_width/2.0,photodet_width/2.0,photodet_thickness/2.0);
    Volume     window_vol(mod_name+"_window", window_box, photodet_mat);

    // update position
    z_placement += photodet_thickness/2.0;
    // place volume
    pv = m_volume.placeVolume(window_vol,Position(0,0,z_placement));
    DetElement   comp_de(mod_de, mod_name + std::string("_sensor_de_") + std::to_string(1) ,  1);
    comp_de.setPlacement(pv);
    // update position
    z_placement += photodet_thickness/2.0;

    // sensitive
    pv.addPhysVolID("sensor", n_sensor);
    window_vol.setSensitiveDetector(sens);
    sensitives[mod_name].push_back(pv);
    ++n_sensor;

    // sensor
    xml_comp_t x_sensor  = x_photodet.child(_Unicode(sensor));
    double     sensor_thickness   = getAttrOrDefault(x_sensor, _U(thickness), 2.0 * mm);
    Material   sensor_mat         = description.material(x_sensor.materialStr());
    int        sensor_nx          = getAttrOrDefault(x_sensor, _Unicode(nx), 2);
    int        sensor_ny          = getAttrOrDefault(x_sensor, _Unicode(ny), 2);

    // layers
    int i_layer = 1;
    for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) {
      xml_comp_t x_layer = li;
      Material layer_mat = description.material(x_layer.materialStr());
      double   layer_thickness = x_layer.thickness();
      Box      layer_box(photodet_width/2.0,photodet_width/2.0,layer_thickness/2.0);
      Volume   layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat);

      // update position
      z_placement += layer_thickness / 2.0;
      // place volume
      pv = m_volume.placeVolume(layer_vol,Position(0,0,z_placement));
      DetElement layer_de(mod_de, mod_name + std::string("_layer_de_") + std::to_string(i_layer),  1);
      layer_de.setPlacement(pv);
      // update position
      z_placement += layer_thickness / 2.0;

      i_layer++;
    }
  }

  //for (size_t ic = 0; ic < sensVols.size(); ++ic) {
  //  PlacedVolume sens_pv = sensVols[ic];
  //  DetElement   comp_de(mod_de, std::string("de_") + sens_pv.volume().name(), ic + 1);
  //  comp_de.setPlacement(sens_pv);
  //  // Acts::ActsExtension* sensorExtension = new Acts::ActsExtension();
  //  //// sensorExtension->addType("sensor", "detector");
  //  // comp_de.addExtension<Acts::ActsExtension>(sensorExtension);
  //  //// comp_de.setAttributes(description, sens_pv.volume(), x_layer.regionStr(),
  //  //// x_layer.limitsStr(),
  //  ////                      xml_det_t(xmleles[m_nam]).visStr());
  //}
  //DetElement window_de(sdet, mod_name + std::string("_window_de") + std::to_string(1), 1);
  //window_de.setPlacement(pv);

  modules[mod_name]                   = m_volume;
  module_assembly_delements[mod_name] = mod_de;
  // end module

  // place modules in the sectors (disk)
  auto points = athena::geo::fillSquares({0., 0.}, mod_width, rmin, rmax);

  // mod_name = ...
  Placements& sensVols = sensitives[mod_name];
  auto        mod_v    = modules[mod_name];
  // determine module direction, always facing z = 0
  double roty = dims.zmin() < 0. ? -M_PI : 0 ;

  // read module positions
  std::vector<std::tuple<double,double,double>> positions;
  for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) {
    xml_comp_t x_positions = x_positions_i;
    for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) {
      xml_comp_t x_position = x_position_i;
      positions.push_back(
        std::make_tuple(x_positions.scale() * x_position.x() * mm,
                        x_positions.scale() * x_position.y() * mm,
                        -x_positions.z0()));
    }
  }
  // if no positions, then autoplacement
  if (positions.empty()) {
    for (double x = mod_width / 2.0; x < rmax - mod_width / 2.0; x += mod_width) {
      for (double y = mod_width / 2.0; y < rmax - mod_width / 2.0; y += mod_width) {
        if (pow(x + mod_width / 2.0,2) + pow(y + mod_width / 2.0,2) > rmax*rmax) continue;
        if (pow(x - mod_width / 2.0,2) + pow(y - mod_width / 2.0,2) < rmin*rmin) continue;
        positions.push_back(std::make_tuple(x, y, 0));
      }
    }
  }

  // place modules
  int i_mod = 1; // starts at 1
  for (auto& p: positions) {

    // get positions in one quadrant
    double x = std::get<0>(p);
    double y = std::get<1>(p);
    double z0 = std::get<2>(p);

    // and place in all quadrants (intentional shadowing)
    for (auto& p: decltype(positions){{x,y,z0}, {y,-x,z0}, {-x,-y,z0}, {-y,x,z0}}) {

      // get positions (intentional shadowing)
      double x = std::get<0>(p);
      double y = std::get<1>(p);
      double z0 = std::get<2>(p);

      Transform3D tr;
      if(projective) {
        double rotAngX = atan(y/z0);
        double rotAngY = -1.*atan(x/z0);
        tr = Translation3D(x, y, 0) * RotationX(rotAngX) * RotationY(rotAngY);
      } else {
        tr = Translation3D(x, y, 0) * RotationX(0);
      }

      // mod placement
      pv = envVol.placeVolume(mod_v, tr);
      pv.addPhysVolID("module", i_mod);

      auto mod_det_element  = module_assembly_delements[mod_name].clone(mod_name + "__" + std::to_string(i_mod));
      mod_det_element.setPlacement(pv);
      sdet.add(mod_det_element);

      i_mod++;
    }
  }

  // additional layers
  if (x_det.hasChild(_Unicode(layer))) {
    xml_comp_t x_layer = x_det.child(_Unicode(layer));
    double   layer_thickness = x_layer.thickness();
    Material layer_mat = description.material(x_layer.materialStr());
    Tube   frameShape(rmin, rmax, layer_thickness / 2., 0., 2 * M_PI);
    Volume frameVol("MRICH_Frame", frameShape, layer_mat);
    pv = envVol.placeVolume(frameVol, Position(0, 0, (length - layer_thickness) / 2.0));
  }

  // place envelope
  Volume motherVol = description.pickMotherVolume(sdet);
  if (reflect) {
    pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, M_PI, 0), Position(0, 0, -zpos)));
  } else {
    pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0,    0, 0), Position(0, 0, +zpos)));
  }
  pv.addPhysVolID("system", x_det.id());
  sdet.setPlacement(pv);
  return sdet;
}


//void addModules(Volume &mother, xml::DetElement &detElem, Detector &description, SensitiveDetector &sens)
//{
//    xml::Component dims = detElem.dimensions();
//    xml::Component mods = detElem.child(_Unicode(modules));
//
//    auto rmin = dims.rmin();
//    auto rmax = dims.rmax();
//
//    auto mThick = mods.attr<double>(_Unicode(thickness));
//    auto mWidth = mods.attr<double>(_Unicode(width));
//    auto mGap = mods.attr<double>(_Unicode(gap));
//
//    auto modMat = description.material(mods.materialStr());
//    auto gasMat = description.material("AirOptical");
//
//    // single module
//    Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm);
//    Volume mVol("ce_MRICH_mod_Solid", mShape, modMat);
//
//    // a thin gas layer to detect optical photons
//    Box modShape(mWidth/2., mWidth/2., mThick/2.);
//    Volume modVol("ce_MRICH_mod_Solid_v", modShape, gasMat);
//    // thin gas layer is on top (+z) of the material
//    modVol.placeVolume(mVol, Position(0., 0., -0.1*mm));
//
//    modVol.setVisAttributes(description.visAttributes(mods.visStr()));
//    sens.setType("tracker");
//    modVol.setSensitiveDetector(sens);
//
//    // place modules in the sectors (disk)
//    auto points = ref::utils::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap);
//
//    // determine module direction, always facing z = 0
//    double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.;
//    int imod = 1;
//    for (auto &p : points) {
//        // operations are inversely ordered
//        Transform3D tr = Translation3D(p.x(), p.y(), 0.)        // move to position
//                       * RotationY(roty);                       // facing z = 0.
//        auto modPV = mother.placeVolume(modVol, tr);
//        modPV.addPhysVolID("sector", 1).addPhysVolID("module", imod ++);
//    }
//}

// clang-format off
DECLARE_DETELEMENT(athena_MRICH, createDetector)