epic/src/ZeroDegreeCalorimeterEcalWSciFi_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Jihee Kim, Wouter Deconinck
0003
0004 #include "DD4hep/DetFactoryHelper.h"
0005 #include "DD4hep/OpticalSurfaces.h"
0006 #include "DD4hep/Printout.h"
0007 #include "DDRec/DetectorData.h"
0008 #include "DDRec/Surface.h"
0009 #include <XML/Helper.h>
0010 #include <algorithm>
0011 #include <iostream>
0012 #include <math.h>
0013 #include <tuple>
0014 //////////////////////////////////////////////////
0015 // Far Forward Ion Zero Degree Calorimeter - Ecal
0016 // Reference from ATHENA ScFiCalorimeter_geo.cpp
0017 //////////////////////////////////////////////////
0018
0019 using namespace std;
0020 using namespace dd4hep;
0021
0022 // main
0023 static Ref_t create_detector(Detector& desc, xml_h handle, SensitiveDetector sens) {
0024   xml::DetElement detElem = handle;
0025   std::string detName     = detElem.nameStr();
0026   int detID               = detElem.id();
0027   DetElement det(detName, detID);
0028   sens.setType("calorimeter");
0029   auto dim      = detElem.dimensions();
0030   auto width    = dim.x();
0031   auto length   = dim.z();
0032   xml_dim_t pos = detElem.position();
0033   xml_dim_t rot = detElem.rotation();
0034
0035   // envelope
0036   Box envShape(width * 0.5, width * 0.5, length * 0.5);
0037   Volume env(detName + "_envelope", envShape, desc.material("Air"));
0038   env.setVisAttributes(desc.visAttributes(detElem.visStr()));
0039
0040   // build module
0041   xml_comp_t mod_x = detElem.child(_Unicode(module));
0042   auto sx          = mod_x.attr<double>(_Unicode(sizex));
0043   auto sy          = mod_x.attr<double>(_Unicode(sizey));
0044   auto sz          = mod_x.attr<double>(_Unicode(sizez));
0045
0046   Box modShape(sx / 2., sy / 2., sz / 2.);
0047   auto modMat = desc.material(mod_x.attr<std::string>(_Unicode(material)));
0048   Volume modVol("module_vol", modShape, modMat);
0049   modVol.setVisAttributes(desc.visAttributes(mod_x.visStr()));
0050   // modVol.setSensitiveDetector(sens);
0051
0052   if (mod_x.hasChild(_Unicode(fiber))) {
0053     auto fiber_x  = mod_x.child(_Unicode(fiber));
0054     auto fr       = fiber_x.attr<double>(_Unicode(radius));
0055     auto fsx      = fiber_x.attr<double>(_Unicode(spacex));
0056     auto fsy      = fiber_x.attr<double>(_Unicode(spacey));
0057     auto foff     = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5 * mm);
0058     auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
0059     Tube fiberShape(0., fr, sz / 2. - 1. * mm);
0060     Volume fiberVol("fiber_vol", fiberShape, fiberMat);
0061     fiberVol.setSensitiveDetector(sens);
0062
0063     // Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
0064     // So each fiber is fully contained in a regular hexagon, which are placed as
0065     // the parameters space x and space y are used to add additional spaces between the hexagons
0066     double fside  = 2. / std::sqrt(3.) * fr;
0067     double fdistx = 2. * fside + fsx;
0068     double fdisty = 2. * fr + fsy;
0069
0070     // maximum numbers of the fibers, help narrow the loop range
0071     int nx = int(sx / (2. * fr)) + 1;
0072     int ny = int(sy / (2. * fr)) + 1;
0073
0074     // place the fibers
0075     double y0   = (foff + fside);
0076     int nfibers = 0;
0077     for (int iy = 0; iy < ny; ++iy) {
0078       double y = y0 + fdisty * iy;
0079       // about to touch the boundary
0080       if ((sy - y) < y0) {
0081         break;
0082       }
0083       double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
0084       for (int ix = 0; ix < nx; ++ix) {
0085         double x = x0 + fdistx * ix;
0086         // about to touch the boundary
0087         if ((sx - x) < x0) {
0088           break;
0089         }
0090         auto fiberPV =
0091             modVol.placeVolume(fiberVol, nfibers++, Position{x - sx / 2., y - sy / 2., 0});
0092         fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
0093       }
0094     }
0095     // if no fibers we make the module itself sensitive
0096   } else {
0097     modVol.setSensitiveDetector(sens);
0098   }
0099
0100   // Module Position
0101   double mod_x_pos = 0.0;
0102   double mod_y_pos = 0.0;
0103   double mod_z_pos = 0.0 * mm;
0104   double mgap      = 0.000001 * mm;
0105   int mNTowers     = floor(width / (sx + mgap));
0106   // std::cout << "mNTowers: " << mNTowers << std::endl;
0107
0108   int k = 0;
0109   // Place Modules
0110   for (int j = 0; j < mNTowers; j++) {
0111     if (j == 0)
0112       mod_y_pos = width * 0.5 - (sy + mgap) * 0.5;
0113     else
0114       mod_y_pos -= (sy + mgap);
0115
0116     for (int i = 0; i < mNTowers; i++) {
0117       if (i == 0)
0118         mod_x_pos = width * 0.5 - (sx + mgap) * 0.5;
0119       else
0120         mod_x_pos -= (sx + mgap);
0121
0122       PlacedVolume pv_mod = env.placeVolume(modVol, Position(mod_x_pos, mod_y_pos, mod_z_pos));
0123       pv_mod.addPhysVolID("module", k++);
0124       // std::cout << "j: " << j << "i: " << i << " Position: " << mod_x_pos << " " << mod_y_pos << " " << mod_z_pos <<
0125       // std::endl;
0126     }
0127   }
0128
0129   // detector position and rotation
0130   Volume motherVol = desc.pickMotherVolume(det);
0131   Transform3D tr(RotationZYX(rot.z(), rot.y(), rot.x()), Position(pos.x(), pos.y(), pos.z()));
0132   PlacedVolume envPV = motherVol.placeVolume(env, tr);
0133   envPV.addPhysVolID("system", detID);
0134   det.setPlacement(envPV);
0135   return det;
0136 }
0137
0138 DECLARE_DETELEMENT(ZDCEcalScFiCalorimeter, create_detector)