Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • EIC/detectors/athena
  • zwzhao/athena
  • FernandoTA/athena
  • palspeic/athena
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1463 additions and 649 deletions
src/HybridCalorimeter_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// A general implementation for homogeneous calorimeter
// it supports three types of placements
// 1. Individual module placement with module dimensions and positions
// 2. Array placement with module dimensions and numbers of row and columns
// 3. Disk placement with module dimensions and (Rmin, Rmax), and (Phimin, Phimax)
// 4. Lines placement with module dimensions and (mirrorx, mirrory)
// (NOTE: anchor point is the 0th block of the line instead of line center)
//--------------------------------------------------------------------------
// Author: Chao Peng (ANL)
// Date: 06/09/2021
//==========================================================================
#include "GeometryHelpers.h"
#include "DD4hep/DetFactoryHelper.h"
#include <XML/Helper.h>
#include <iostream>
#include <algorithm>
#include <tuple>
#include <math.h>
#include <fmt/core.h>
using namespace dd4hep;
// main
static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
{
using namespace std;
using namespace fmt;
xml::DetElement detElem = handle;
std::string detName = detElem.nameStr();
int detID = detElem.id();
DetElement det(detName, detID);
sens.setType("calorimeter");
auto glass_material = desc.material("SciGlass");
auto crystal_material = desc.material("PbWO4");
auto air_material = desc.material("Air");
double ROut = desc.constantAsDouble("EcalEndcapN_rmax");
double RIn_el = desc.constantAsDouble("EcalEndcapN_rmin1");
double ionCutout_dphi = desc.constantAsDouble("EcalEndcapNIonCutout_dphi");
double RIn = desc.constantAsDouble("EcalEndcapN_rmin2");
double SizeZ = desc.constantAsDouble("EcalEndcapN_thickness");
double thickness = desc.constantAsDouble("EcalEndcapN_thickness");
double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
double Glass_z0 = desc.constantAsDouble("GlassModule_z0");
double Glass_Width = desc.constantAsDouble("GlassModule_width");
double Glass_thickness = desc.constantAsDouble("GlassModule_length");
double Glass_Gap = desc.constantAsDouble("GlassModule_wrap");
double glass_distance = desc.constantAsDouble("GlassModule_distance");
double Crystal_Width = desc.constantAsDouble("CrystalModule_width");
double Crystal_thickness = desc.constantAsDouble("CrystalModule_length");
double Crystal_Gap = desc.constantAsDouble("CrystalModule_wrap");
double crystal_distance = desc.constantAsDouble("CrystalModule_distance");
double Crystal_z0 = desc.constantAsDouble("CrystalModule_z0");
// RIn and ROut will define outer tube embedding the calorimeter
// centers_rmin/out define the maximum radius of module centers
// so that modules are not overlapping with mother tube volume
const double glassHypotenuse = std::hypot(glass_distance, glass_distance)/2;
const double crystalHypotenuse = glassHypotenuse/2;
// Offset these values a bit so we don't miss edge-blocks
const double glassCenters_rmax = ROut - glassHypotenuse + 1 * mm;
const double crystalCenters_rmin = RIn + crystalHypotenuse - 1 * mm;
// Also limits of the inner crystal blocks fill
const double cutout_tan = tan(ionCutout_dphi/2);
const double cutout_rmin = RIn_el + crystalHypotenuse - 1 * mm;
// Offset to align the modules at the zmin of the endcap,
const double Crystal_offset = -0.5 * (Crystal_thickness - thickness);
const double Glass_offset = -0.5*(Glass_thickness - thickness);
// envelope
// consists of an glass tube of the full thickness, and a crystal inner tube
// for the crystal that allows us to get closet to the beampipe without
// overlaps, and a partial electron tube that allows us to get closer to the
// electron beampipe in the region where there is no ion beampipe
Tube glass_tube(min(RIn + glassHypotenuse*2, trans_radius), ROut, SizeZ / 2.0, 0., 360.0 * deg);
Tube crystal_tube(RIn, min(RIn + glassHypotenuse*2, trans_radius), Crystal_thickness/ 2.0, 0., 360.0 * deg);
Tube electron_tube(RIn_el, RIn, Crystal_thickness / 2., ionCutout_dphi / 2., 360.0 * deg - ionCutout_dphi / 2);
UnionSolid outer_envelope(glass_tube, crystal_tube, Position(0,0,Crystal_offset));
UnionSolid envelope(outer_envelope, electron_tube, Position(0,0,Crystal_offset));
Volume ecal_vol("negative_ecal", envelope, air_material);
ecal_vol.setVisAttributes(desc.visAttributes("HybridEcalOuterVis"));
// TODO why 1cm and not something else?
double Glass_OuterR = ROut - 1 * cm ;
double Glass_InnerR = trans_radius;
// Geometry of modules
Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_thickness * 0.5);
Volume glass_module("glass_module", glass_box, glass_material);
glass_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
glass_module.setSensitiveDetector(sens);
Box crystal_box("crystal_box", Crystal_Width* 0.5, Crystal_Width * 0.5, Crystal_thickness * 0.5);
Volume crystal_module("crystal_module", crystal_box, crystal_material);
crystal_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
crystal_module.setSensitiveDetector(sens);
// GLASS
double diameter = 2 * Glass_OuterR;
// Can we fit an even or odd amount of glass blocks within our rmax?
// This determines the transition points between crystal and glass as we need the
// outer crystal arangement to be in multiples of 2 (aligned with glass)
const int towersInRow = std::ceil((diameter + Glass_Gap) / (Glass_Width + Glass_Gap));
const bool align_even = (towersInRow % 2);
// Is it odd or even number of towersInRow
double leftTowerPos, topTowerPos;
if(towersInRow%2) {
// |
// [ ][ ][ ][ ][ ]
// ^ |
int towersInHalfRow = std::ceil(towersInRow/2.0);
topTowerPos = leftTowerPos = -towersInHalfRow * (Glass_Width + Glass_Gap);
} else {
// |
// [ ][ ][ ][ ][ ][ ]
// ^ |
int towersInHalfRow = towersInRow/2;
topTowerPos = leftTowerPos = -(towersInHalfRow - 0.5) * (Glass_Width + Glass_Gap);
}
int moduleIndex = 0;
int glass_module_index = 0;
int cryst_module_index = 0;
for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
for(int colIndex=0; colIndex < towersInRow; colIndex++) {
const double glass_x = leftTowerPos + colIndex * (Glass_Width + Glass_Gap);
const double glass_y = topTowerPos + rowIndex * (Glass_Width + Glass_Gap);
const double r = std::hypot(glass_x, glass_y);
// crystal if within the transition radius (as defined by the equivalent glass
// block)
if (r < trans_radius) {
for (const auto dx : {-1, 1}) {
for (const auto& dy : {-1, 1}) {
const double crystal_x = glass_x + dx * crystal_distance / 2;
const double crystal_y = glass_y + dy * crystal_distance / 2;
const double crystal_r = std::hypot(crystal_x, crystal_y);
// check if crystal in the main crystal ring?
const bool mainRing = (crystal_r > crystalCenters_rmin);
const bool innerRing = !mainRing && crystal_r > cutout_rmin;
const bool ionCutout = crystal_x > 0 && fabs(crystal_y / crystal_x) < cutout_tan;
if (mainRing || (innerRing && !ionCutout)) {
auto placement =
ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
placement.addPhysVolID("sector", 1);
placement.addPhysVolID("module", cryst_module_index++);
}
}
}
// Glass block if within the rmax
} else if (r < glassCenters_rmax) {
// glass module
auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, Glass_z0 + Glass_offset));
placement.addPhysVolID("sector", 2);
placement.addPhysVolID("module", glass_module_index++);
}
}
}
desc.add(Constant("EcalEndcapN_NModules_Sector1", std::to_string(cryst_module_index)));
desc.add(Constant("EcalEndcapN_NModules_Sector2", std::to_string(glass_module_index)));
// fmt::print("Total Glass modules: {}\n", towerIndex);
// fmt::print("CE EMCAL GLASS END\n\n");
// detector position and rotation
auto pos = detElem.position();
auto rot = detElem.rotation();
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(rot.z(), rot.y(), rot.x()), Position(pos.x(), pos.y(), pos.z()));
PlacedVolume envPV = motherVol.placeVolume(ecal_vol, tr);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
//@}
DECLARE_DETELEMENT(HybridCalorimeter, create_detector)
src/MRich_geo.cpp 0 → 100644
View file @ 1008a184
//
// 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)
src/PFRICH_geo.cpp 0 → 100644
View file @ 1008a184
//----------------------------------
// pfRICH: Proximity Focusing RICH
// Author: C. Dilks
//----------------------------------
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/Printout.h"
#include "DDRec/DetectorData.h"
#include "DDRec/Surface.h"
#include "GeometryHelpers.h"
#include "Math/Point2D.h"
#include "TMath.h"
#include "TString.h"
#include <XML/Helper.h>
using namespace dd4hep;
using namespace dd4hep::rec;
// 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();
OpticalSurfaceManager surfMgr = desc.surfaceManager();
// attributes -----------------------------------------------------------
// - vessel
double vesselLength = dims.attr<double>(_Unicode(length));
double vesselZmin = dims.attr<double>(_Unicode(zmin));
double vesselZmax = dims.attr<double>(_Unicode(zmax));
double vesselRmin0 = dims.attr<double>(_Unicode(rmin0));
double vesselRmin1 = dims.attr<double>(_Unicode(rmin1));
double vesselRmax0 = dims.attr<double>(_Unicode(rmax0));
double vesselRmax1 = dims.attr<double>(_Unicode(rmax1));
double wallThickness = dims.attr<double>(_Unicode(wall_thickness));
double windowThickness = dims.attr<double>(_Unicode(window_thickness));
auto vesselMat = desc.material(detElem.attr<std::string>(_Unicode(material)));
auto gasvolMat = desc.material(detElem.attr<std::string>(_Unicode(gas)));
auto vesselVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
auto gasvolVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
// - radiator (applies to aerogel and filter)
auto radiatorElem = detElem.child(_Unicode(radiator));
double radiatorRmin = radiatorElem.attr<double>(_Unicode(rmin));
double radiatorRmax = radiatorElem.attr<double>(_Unicode(rmax));
double radiatorPhiw = radiatorElem.attr<double>(_Unicode(phiw));
double radiatorPitch = radiatorElem.attr<double>(_Unicode(pitch));
double radiatorFrontplane = radiatorElem.attr<double>(_Unicode(frontplane));
// - aerogel
auto aerogelElem = radiatorElem.child(_Unicode(aerogel));
auto aerogelMat = desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
auto aerogelVis = desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
double aerogelThickness = aerogelElem.attr<double>(_Unicode(thickness));
// - filter
auto filterElem = radiatorElem.child(_Unicode(filter));
auto filterMat = desc.material(filterElem.attr<std::string>(_Unicode(material)));
auto filterVis = desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
double filterThickness = filterElem.attr<double>(_Unicode(thickness));
// - sensor module
auto sensorElem = detElem.child(_Unicode(sensors)).child(_Unicode(module));
auto sensorMat = desc.material(sensorElem.attr<std::string>(_Unicode(material)));
auto sensorVis = desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
auto sensorSurf = surfMgr.opticalSurface(sensorElem.attr<std::string>(_Unicode(surface)));
double sensorSide = sensorElem.attr<double>(_Unicode(side));
double sensorGap = sensorElem.attr<double>(_Unicode(gap));
double sensorThickness = sensorElem.attr<double>(_Unicode(thickness));
// - sensor plane
auto sensorPlaneElem = detElem.child(_Unicode(sensors)).child(_Unicode(plane));
double sensorPlaneDist = sensorPlaneElem.attr<double>(_Unicode(sensordist));
double sensorPlaneRmin = sensorPlaneElem.attr<double>(_Unicode(rmin));
double sensorPlaneRmax = sensorPlaneElem.attr<double>(_Unicode(rmax));
// - debugging switches
int debug_optics_mode = detElem.attr<int>(_Unicode(debug_optics));
// if debugging optics, override some settings
bool debug_optics = debug_optics_mode > 0;
if (debug_optics) {
printout(WARNING, "PFRICH_geo", "DEBUGGING PFRICH OPTICS");
switch (debug_optics_mode) {
case 1:
vesselMat = aerogelMat = filterMat = sensorMat = gasvolMat = desc.material("VacuumOptical");
break;
case 2:
vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical");
break;
default:
printout(FATAL, "PFRICH_geo", "UNKNOWN debug_optics_mode");
return det;
};
aerogelVis = sensorVis;
gasvolVis = vesselVis = desc.invisible();
};
// BUILD VESSEL //////////////////////////////////////
/* - `vessel`: aluminum enclosure, the mother volume of the pfRICH
* - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
* are children of `gasvol`
*/
// tank solids
double boreDelta = vesselRmin1 - vesselRmin0;
Cone vesselTank(vesselLength / 2.0, vesselRmin1, vesselRmax1, vesselRmin0, vesselRmax0);
Cone gasvolTank(vesselLength / 2.0 - windowThickness, vesselRmin1 + wallThickness, vesselRmax1 - wallThickness,
vesselRmin0 + wallThickness, vesselRmax0 - wallThickness);
// extra solids for `debug_optics` only
Box vesselBox(1001, 1001, 1001);
Box gasvolBox(1000, 1000, 1000);
// choose vessel and gasvol solids (depending on `debug_optics_mode` (0=disabled))
Solid vesselSolid, gasvolSolid;
switch (debug_optics_mode) {
case 0:
vesselSolid = vesselTank;
gasvolSolid = gasvolTank;
break; // `!debug_optics`
case 1:
vesselSolid = vesselBox;
gasvolSolid = gasvolBox;
break;
case 2:
vesselSolid = vesselBox;
gasvolSolid = gasvolTank;
break;
};
// volumes
Volume vesselVol(detName, vesselSolid, vesselMat);
Volume gasvolVol(detName + "_gas", gasvolSolid, gasvolMat);
vesselVol.setVisAttributes(vesselVis);
gasvolVol.setVisAttributes(gasvolVis);
// reference positions
// - the vessel is created such that the center of the cylindrical tank volume
// coincides with the origin; this is called the "origin position" of the vessel
// - when the vessel (and its children volumes) is placed, it is translated in
// the z-direction to be in the proper ATHENA-integration location
// - these reference positions are for the frontplane and backplane of the vessel,
// with respect to the vessel origin position
auto originFront = Position(0., 0., vesselLength / 2.0);
auto originBack = Position(0., 0., -vesselLength / 2.0);
// sensitive detector type
sens.setType("tracker");
// BUILD RADIATOR //////////////////////////////////////
// attributes
double airGap = 0.01*mm; // air gap between aerogel and filter (FIXME? actually it's currently a gas gap)
// solid and volume: create aerogel and filter
Cone aerogelSolid(
aerogelThickness/2,
radiatorRmin + boreDelta * aerogelThickness / vesselLength, /* at backplane */
radiatorRmax,
radiatorRmin, /* at frontplane */
radiatorRmax
);
Cone filterSolid(
filterThickness/2,
radiatorRmin + boreDelta * (aerogelThickness + airGap + filterThickness) / vesselLength, /* at backplane */
radiatorRmax,
radiatorRmin + boreDelta * (aerogelThickness + airGap) / vesselLength, /* at frontplane */
radiatorRmax
);
Volume aerogelVol(detName + "_aerogel", aerogelSolid, aerogelMat);
Volume filterVol(detName + "_filter", filterSolid, filterMat);
aerogelVol.setVisAttributes(aerogelVis);
filterVol.setVisAttributes(filterVis);
// aerogel placement and surface properties
// TODO [low-priority]: define skin properties for aerogel and filter
auto radiatorPos = Position(0., 0., radiatorFrontplane - 0.5 * aerogelThickness) + originFront;
auto aerogelPV = gasvolVol.placeVolume(
aerogelVol,
Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
* RotationY(radiatorPitch) // change polar angle to specified pitch // (FIXME: probably broken, currently not in use)
);
DetElement aerogelDE(det, "aerogel_de", 0);
aerogelDE.setPlacement(aerogelPV);
// SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
// aerogelSkin.isValid();
// filter placement and surface properties
if (!debug_optics) {
auto filterPV = gasvolVol.placeVolume(
filterVol,
Translation3D(0., 0., -airGap) // add an airgap (FIXME: actually a gas gap)
* Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
* RotationY(radiatorPitch) // change polar angle
* Translation3D(0., 0., -(aerogelThickness + filterThickness) / 2.) // move to aerogel backplane
);
DetElement filterDE(det, "filter_de", 0);
filterDE.setPlacement(filterPV);
// SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
// filterSkin.isValid();
};
// BUILD SENSORS ///////////////////////
// solid and volume: single sensor module
Box sensorSolid(sensorSide / 2., sensorSide / 2., sensorThickness / 2.);
Volume sensorVol(detName + "_sensor", sensorSolid, sensorMat);
sensorVol.setVisAttributes(sensorVis);
// sensitivity
if (!debug_optics)
sensorVol.setSensitiveDetector(sens);
// sensor plane positioning: we want `sensorPlaneDist` to be the distance between the
// aerogel backplane (i.e., aerogel/filter boundary) and the sensor active surface (e.g, photocathode)
double sensorZpos = radiatorFrontplane - aerogelThickness - sensorPlaneDist - 0.5 * sensorThickness;
auto sensorPlanePos = Position(0., 0., sensorZpos) + originFront; // reference position
// miscellaneous
int imod = 0; // module number
double tBoxMax = vesselRmax1; // sensors will be tiled in tBox, within annular limits
// SENSOR MODULE LOOP ------------------------
/* cartesian tiling loop
* - start at (x=0,y=0), to center the grid
* - loop over positive-x positions; for each, place the corresponding negative-x sensor too
* - nested similar loop over y positions
*/
double sx, sy;
for (double usx = 0; usx <= tBoxMax; usx += sensorSide + sensorGap) {
for (int sgnx = 1; sgnx >= (usx > 0 ? -1 : 1); sgnx -= 2) {
for (double usy = 0; usy <= tBoxMax; usy += sensorSide + sensorGap) {
for (int sgny = 1; sgny >= (usy > 0 ? -1 : 1); sgny -= 2) {
// sensor (x,y) center
sx = sgnx * usx;
sy = sgny * usy;
// annular cut
if (std::hypot(sx, sy) < sensorPlaneRmin || std::hypot(sx, sy) > sensorPlaneRmax)
continue;
// placement (note: transformations are in reverse order)
auto sensorPV = gasvolVol.placeVolume(
sensorVol, Transform3D(Translation3D(sensorPlanePos.x(), sensorPlanePos.y(),
sensorPlanePos.z()) // move to reference position
* Translation3D(sx, sy, 0.) // move to grid position
));
// generate LUT for module number -> sensor position, for readout mapping tests
// printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
// properties
sensorPV.addPhysVolID("module", imod);
DetElement sensorDE(det, Form("sensor_de_%d", imod), imod);
sensorDE.setPlacement(sensorPV);
if (!debug_optics) {
SkinSurface sensorSkin(desc, sensorDE, "sensor_optical_surface", sensorSurf,
sensorVol); // TODO: 3rd arg needs `imod`?
sensorSkin.isValid();
};
// increment sensor module number
imod++;
};
};
};
};
// END SENSOR MODULE LOOP ------------------------
//
// Add service material if desired
if (detElem.child("sensors").hasChild("services")) {
xml_comp_t x_service = detElem.child("sensors").child(_Unicode(services));
Assembly service_vol("services");
service_vol.setVisAttributes(desc, x_service.visStr());
// Compute service total thickness from components
double total_thickness = 0;
xml_coll_t ci(x_service, _Unicode(component));
for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
total_thickness += xml_comp_t(ci).thickness();
}
int ncomponents = 0;
double thickness_sum = -total_thickness / 2.0;
for (xml_coll_t ci(x_service, _Unicode(component)); ci; ++ci, ncomponents++) {
xml_comp_t x_comp = ci;
double thickness = x_comp.thickness();
Tube c_tube{sensorPlaneRmin, sensorPlaneRmax, thickness/2};
Volume c_vol{_toString(ncomponents, "component%d"), c_tube, desc.material(x_comp.materialStr())};
c_vol.setVisAttributes(desc, x_comp.visStr());
service_vol.placeVolume(c_vol, Position(0, 0, thickness_sum + thickness / 2.0));
thickness_sum += thickness;
}
gasvolVol.placeVolume(service_vol,
Transform3D(Translation3D(sensorPlanePos.x(), sensorPlanePos.y(),
sensorPlanePos.z() - sensorThickness - total_thickness)));
}
// place gas volume
PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol, Position(0, 0, 0));
DetElement gasvolDE(det, "gasvol_de", 0);
gasvolDE.setPlacement(gasvolPV);
// place mother volume (vessel)
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume vesselPV = motherVol.placeVolume(vesselVol, Position(0, 0, vesselZmin) - originFront);
vesselPV.addPhysVolID("system", detID);
det.setPlacement(vesselPV);
return det;
};
// clang-format off
DECLARE_DETELEMENT(athena_PFRICH, createDetector)
src/PolyhedraEndcapCalorimeter2_geo.cpp
View file @ 1008a184
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
// //
//========================================================================== //==========================================================================
// //
// Modified for TOPSiDE detector // Modified for ATHENA detector
// //
//========================================================================== //==========================================================================
#include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DetFactoryHelper.h"
...@@ -32,7 +32,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -32,7 +32,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
int numsides = dim.numsides(); int numsides = dim.numsides();
xml::Component pos = x_det.position(); xml::Component pos = x_det.position();
double rmin = dim.rmin(); double rmin = dim.rmin();
double rmax = dim.rmax() * std::cos(M_PI / numsides); double rmax = dim.rmax();
double zmin = dim.zmin(); double zmin = dim.zmin();
Layering layering(x_det); Layering layering(x_det);
double totalThickness = layering.totalThickness(); double totalThickness = layering.totalThickness();
...@@ -128,4 +128,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -128,4 +128,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
} }
// clang-format off // clang-format off
DECLARE_DETELEMENT(refdet_PolyhedraEndcapCalorimeter2, create_detector) DECLARE_DETELEMENT(athena_PolyhedraEndcapCalorimeter2, create_detector)
DECLARE_DETELEMENT(athena_PolyhedraEndcapCalorimeter, create_detector)
src/ScFiCalorimeter_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// Scintillating fiber calorimeter with tower shape blocks
// reference: https://github.com/adamjaro/lmon/blob/master/calo/src/WScFiZXv3.cxx
// Support disk placement
//--------------------------------------------------------------------------
// Author: Chao Peng (ANL)
// Date: 07/19/2021
//==========================================================================
#include "GeometryHelpers.h"
#include "DD4hep/DetFactoryHelper.h"
#include <XML/Helper.h>
#include <iostream>
#include <algorithm>
#include <tuple>
#include <math.h>
using namespace dd4hep;
using Point = ROOT::Math::XYPoint;
std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens);
// helper function to get x, y, z if defined in a xml component
template<class XmlComp>
Position get_xml_xyz(const XmlComp &comp, dd4hep::xml::Strng_t name)
{
Position pos(0., 0., 0.);
if (comp.hasChild(name)) {
auto child = comp.child(name);
pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
}
return pos;
}
// main
static Ref_t create_detector(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);
sens.setType("calorimeter");
auto dim = detElem.dimensions();
auto rmin = dim.rmin();
auto rmax = dim.rmax();
auto length = dim.length();
auto phimin = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimin), 0.);
auto phimax = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimax), 2.*M_PI);
// envelope
Tube envShape(rmin, rmax, length/2., phimin, phimax);
Volume env(detName + "_envelope", envShape, desc.material("Air"));
env.setVisAttributes(desc.visAttributes(detElem.visStr()));
// build module
auto [modVol, modSize] = build_module(desc, detElem.child(_Unicode(module)), sens);
double modSizeR = std::sqrt(modSize.x() * modSize.x() + modSize.y() * modSize.y());
double assembly_rwidth = modSizeR*2.;
int nas = int((rmax - rmin) / assembly_rwidth) + 1;
std::vector<Assembly> assemblies;
// calorimeter block z-offsets (as blocks are shorter than the volume length)
const double block_offset = -0.5*(length - modSize.z());
for (int i = 0; i < nas; ++i) {
Assembly assembly(detName + Form("_ring%d", i + 1));
auto assemblyPV = env.placeVolume(assembly, Position{0., 0., block_offset});
assemblyPV.addPhysVolID("ring", i + 1);
assemblies.emplace_back(std::move(assembly));
}
// std::cout << assemblies.size() << std::endl;
int modid = 1;
for (int ix = 0; ix < int(2.*rmax / modSize.x()) + 1; ++ix) {
double mx = modSize.x() * ix - rmax;
for (int iy = 0; iy < int(2.*rmax / modSize.y()) + 1; ++iy) {
double my = modSize.y() * iy - rmax;
double mr = std::sqrt(mx*mx + my*my);
if (mr - modSizeR >= rmin && mr + modSizeR <= rmax) {
int ias = int((mr - rmin) / assembly_rwidth);
auto &assembly = assemblies[ias];
auto modPV = assembly.placeVolume(modVol, Position(mx, my, 0.));
modPV.addPhysVolID("module", modid++);
}
}
}
desc.add(Constant(detName + "_NModules", std::to_string(modid - 1)));
for (auto &assembly : assemblies) {
assembly.ptr()->Voxelize("");
}
// detector position and rotation
auto pos = get_xml_xyz(detElem, _Unicode(position));
auto rot = get_xml_xyz(detElem, _Unicode(rotation));
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
PlacedVolume envPV = motherVol.placeVolume(env, tr);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// helper function to build module with scintillating fibers
std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens)
{
auto sx = mod_x.attr<double>(_Unicode(sizex));
auto sy = mod_x.attr<double>(_Unicode(sizey));
auto sz = mod_x.attr<double>(_Unicode(sizez));
Box modShape(sx/2., sy/2., sz/2.);
auto modMat = desc.material(mod_x.attr<std::string>(_Unicode(material)));
Volume modVol("module_vol", modShape, modMat);
if (mod_x.hasAttr(_Unicode(vis))) {
modVol.setVisAttributes(desc.visAttributes(mod_x.attr<std::string>(_Unicode(vis))));
}
if (mod_x.hasChild("fiber")) {
auto fiber_x = mod_x.child(_Unicode(fiber));
auto fr = fiber_x.attr<double>(_Unicode(radius));
auto fsx = fiber_x.attr<double>(_Unicode(spacex));
auto fsy = fiber_x.attr<double>(_Unicode(spacey));
auto foff = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5*mm);
auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
Tube fiberShape(0., fr, sz/2.);
Volume fiberVol("fiber_vol", fiberShape, fiberMat);
fiberVol.setSensitiveDetector(sens);
// Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
// So each fiber is fully contained in a regular hexagon, which are placed as
// ______________________________________
// | ____ ____ |
// even: | / \ / \ |
// | ____/ \____/ \____ |
// | / \ / \ / \ |
// odd: | / \____/ \____/ \ |
// | \ / \ / \ / |
// | \____/ \____/ \____/ |
// even: | \ / \ / |
// | \____/ \____/ ___|___
// |____________________________________|___offset
// | |
// |offset
// the parameters space x and space y are used to add additional spaces between the hexagons
double fside = 2. / std::sqrt(3.) * fr;
double fdistx = 2. * fside + fsx;
double fdisty = 2. * fr + fsy;
// maximum numbers of the fibers, help narrow the loop range
int nx = int(sx / (2.*fr)) + 1;
int ny = int(sy / (2.*fr)) + 1;
// std::cout << sx << ", " << sy << ", " << fr << ", " << nx << ", " << ny << std::endl;
// place the fibers
double y0 = (foff + fside);
int nfibers = 0;
for (int iy = 0; iy < ny; ++iy) {
double y = y0 + fdisty * iy;
// about to touch the boundary
if ((sy - y) < y0) { break; }
double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
for (int ix = 0; ix < nx; ++ix) {
double x = x0 + fdistx * ix;
// about to touch the boundary
if ((sx - x) < x0) { break; }
auto fiberPV = modVol.placeVolume(fiberVol, nfibers++, Position{x - sx/2., y - sy/2., 0});
//std::cout << "(" << ix << ", " << iy << ", " << x - sx/2. << ", " << y - sy/2. << ", " << fr << "),\n";
fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
}
}
// if no fibers we make the module itself sensitive
} else {
modVol.setSensitiveDetector(sens);
}
return std::make_tuple(modVol, Position{sx, sy, sz});
}
DECLARE_DETELEMENT(ScFiCalorimeter, create_detector)
src/SimpleDiskTracker_geo.cpp → src/SimpleDiskDetector_geo.cpp
View file @ 1008a184
...@@ -15,14 +15,20 @@ ...@@ -15,14 +15,20 @@
// //
//========================================================================== //==========================================================================
#include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DetFactoryHelper.h"
#if defined(USE_ACTSDD4HEP)
#include "ActsDD4hep/ActsExtension.hpp"
#include "ActsDD4hep/ConvertMaterial.hpp"
#else
#include "Acts/Plugins/DD4hep/ActsExtension.hpp" #include "Acts/Plugins/DD4hep/ActsExtension.hpp"
#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp" #include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
#endif
using namespace std; using namespace std;
using namespace dd4hep; using namespace dd4hep;
using namespace dd4hep::detail; using namespace dd4hep::detail;
static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) static Ref_t SimpleDiskDetector_create_detector(Detector& description, xml_h e, SensitiveDetector sens)
{ {
xml_det_t x_det = e; xml_det_t x_det = e;
Material air = description.air(); Material air = description.air();
...@@ -126,6 +132,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -126,6 +132,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
return sdet; return sdet;
} }
DECLARE_DETELEMENT(athena_SimpleDiskTracker, create_detector) DECLARE_DETELEMENT(ref_SolenoidEndcap, SimpleDiskDetector_create_detector)
DECLARE_DETELEMENT(ref_DiskTracker, create_detector) DECLARE_DETELEMENT(athena_SolenoidEndcap, SimpleDiskDetector_create_detector)
DECLARE_DETELEMENT(ref_SolenoidEndcap, create_detector)
src/SimpleRectangularTracker_geo.cpp deleted 100644 → 0
View file @ 2e7fa037
//==========================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author : M.Frank
// Modified : W.Armstrong
//
//==========================================================================
//
// Specialized generic detector constructor
//
//==========================================================================
#include "DD4hep/DetFactoryHelper.h"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::detail;
static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
xml_det_t x_det = e;
Material air = description.air();
string det_name = x_det.nameStr();
bool reflect = x_det.reflect();
DetElement sdet(det_name,x_det.id());
Assembly assembly(det_name);
PlacedVolume pv;
int l_num = 0;
xml::Component pos = x_det.position();
for(xml_coll_t i(x_det,_U(layer)); i; ++i, ++l_num) {
xml_comp_t x_layer = i;
string l_nam = det_name + _toString(l_num, "_layer%d");
double x_lay = x_layer.x();
double y_lay = x_layer.y();
double z = 0;
double zmin = 0;
double layerWidth = 0.;
int s_num = 0;
for(xml_coll_t j(x_layer,_U(slice)); j; ++j) {
double thickness = xml_comp_t(j).thickness();
layerWidth += thickness;
}
Box l_box(x_lay/2.0, y_lay/2.0, layerWidth/2.0 );
Volume l_vol(l_nam, l_box, air);
l_vol.setVisAttributes(description, x_layer.visStr());
for (xml_coll_t j(x_layer, _U(slice)); j; ++j, ++s_num) {
xml_comp_t x_slice = j;
double thick = x_slice.thickness();
Material mat = description.material(x_slice.materialStr());
string s_nam = l_nam + _toString(s_num, "_slice%d");
Volume s_vol(s_nam, Box(x_lay/2.0, y_lay/2.0, thick/2.0), mat);
if (x_slice.isSensitive()) {
sens.setType("tracker");
s_vol.setSensitiveDetector(sens);
}
s_vol.setAttributes(description, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
pv = l_vol.placeVolume(s_vol, Position(0, 0, z - zmin - layerWidth / 2 + thick / 2));
pv.addPhysVolID("slice", s_num);
}
DetElement layer(sdet,l_nam+"_pos",l_num);
pv = assembly.placeVolume(l_vol,Position(0,0,zmin+layerWidth/2.));
pv.addPhysVolID("layer",l_num);
pv.addPhysVolID("barrel",1);
layer.setPlacement(pv);
if ( reflect ) {
pv = assembly.placeVolume(l_vol,Transform3D(RotationY(M_PI),Position(0,0,-zmin-layerWidth/2)));
pv.addPhysVolID("layer",l_num);
pv.addPhysVolID("barrel",2);
DetElement layerR = layer.clone(l_nam+"_neg");
sdet.add(layerR.setPlacement(pv));
}
}
if ( x_det.hasAttr(_U(combineHits)) ) {
sdet.setCombineHits(x_det.attr<bool>(_U(combineHits)),sens);
}
pv = description.pickMotherVolume(sdet).placeVolume(assembly,Position(pos.x(),pos.y(),pos.z()));
pv.addPhysVolID("system", x_det.id()); // Set the subdetector system ID.
sdet.setPlacement(pv);
return sdet;
}
DECLARE_DETELEMENT(ref_RectangularTracker,create_detector)
src/ref_SolenoidCoil_geo.cpp → src/SolenoidCoil_geo.cpp
View file @ 1008a184
...@@ -89,4 +89,4 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -89,4 +89,4 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
return sdet; return sdet;
} }
DECLARE_DETELEMENT(refdet_SolenoidCoil,create_detector) DECLARE_DETELEMENT(athena_SolenoidCoil,create_detector)
src/TrapEndcapTracker_geo.cpp
View file @ 1008a184
//========================================================================== /** \addtogroup Trackers Trackers
// AIDA Detector description implementation * \brief Type: **BarrelTrackerWithFrame**.
//-------------------------------------------------------------------------- * \author W. Armstrong
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) *
// All rights reserved. * \ingroup trackers
// *
// For the licensing terms see $DD4hepINSTALL/LICENSE. * @{
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. */
// #include <array>
// Author : M.Frank
//
//==========================================================================
//
// Specialized generic detector constructor
//
//==========================================================================
#include <map> #include <map>
#include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
#include "DD4hep/Shapes.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "XML/Utilities.h"
#include "XML/Layering.h"
#if defined(USE_ACTSDD4HEP)
#include "ActsDD4hep/ActsExtension.hpp"
#include "ActsDD4hep/ConvertMaterial.hpp"
#else
#include "Acts/Plugins/DD4hep/ActsExtension.hpp" #include "Acts/Plugins/DD4hep/ActsExtension.hpp"
#include "Acts/Definitions/Units.hpp" #include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
#endif
using namespace std; using namespace std;
using namespace dd4hep; using namespace dd4hep;
using namespace dd4hep::rec;
using namespace dd4hep::detail; using namespace dd4hep::detail;
/** Endcap Trapezoidal Tracker.
*
* @author Whitney Armstrong
*
*/
static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)
{ {
typedef vector<PlacedVolume> Placements; typedef vector<PlacedVolume> Placements;
...@@ -34,23 +45,74 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -34,23 +45,74 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
DetElement sdet(det_name, det_id); DetElement sdet(det_name, det_id);
Assembly assembly(det_name); Assembly assembly(det_name);
Material air = description.material("Air"); Material air = description.material("Air");
// Volume assembly (det_name,Box(10000,10000,10000),vacuum); // Volume assembly (det_name,Box(10000,10000,10000),vacuum);
Volume motherVol = description.pickMotherVolume(sdet); Volume motherVol = description.pickMotherVolume(sdet);
int m_id = 0, c_id = 0, n_sensor = 0; int m_id = 0, c_id = 0, n_sensor = 0;
map<string, Volume> modules; map<string, Volume> modules;
map<string, Placements> sensitives; map<string, Placements> sensitives;
map<string, std::vector<VolPlane>> volplane_surfaces;
map<string, std::array<double, 2>> module_thicknesses;
PlacedVolume pv; PlacedVolume pv;
Acts::ActsExtension* detWorldExt = new Acts::ActsExtension(); // ACTS extension
detWorldExt->addType("endcap", "detector"); {
sdet.addExtension<Acts::ActsExtension>(detWorldExt); Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
detWorldExt->addType("endcap", "detector");
// SJJ probably need to set the envelope here, as ACTS can't figure
// that out for Assembly volumes. May also need binning to properly pick up
// on the support material @TODO
//
// Add the volume boundary material if configured
for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
xml_comp_t x_boundary_material = bmat;
Acts::xmlToProtoSurfaceMaterial(x_boundary_material, *detWorldExt, "boundary_material");
}
sdet.addExtension<Acts::ActsExtension>(detWorldExt);
}
assembly.setVisAttributes(description.invisible()); assembly.setVisAttributes(description.invisible());
sens.setType("tracker"); sens.setType("tracker");
for (xml_coll_t su(x_det, _U(support)); su; ++su) {
xml_comp_t x_support = su;
double support_thickness = getAttrOrDefault(x_support, _U(thickness), 2.0 * mm);
double support_length = getAttrOrDefault(x_support, _U(length), 2.0 * mm);
double support_rmin = getAttrOrDefault(x_support, _U(rmin), 2.0 * mm);
double support_zstart = getAttrOrDefault(x_support, _U(zstart), 2.0 * mm);
std::string support_name = getAttrOrDefault<std::string>(x_support, _Unicode(name), "support_tube");
std::string support_vis = getAttrOrDefault<std::string>(x_support, _Unicode(vis), "AnlRed");
xml_dim_t pos (x_support.child(_U(position), false));
xml_dim_t rot (x_support.child(_U(rotation), false));
Solid support_solid;
if(x_support.hasChild("shape")){
xml_comp_t shape(x_support.child(_U(shape)));
string shape_type = shape.typeStr();
support_solid = xml::createShape(description, shape_type, shape);
} else {
support_solid = Tube(support_rmin, support_rmin + support_thickness, support_length / 2);
}
Transform3D tr = Transform3D(Rotation3D(),Position(0,0,(reflect?-1.0:1.0) * (support_zstart + support_length / 2)));
if ( pos.ptr() && rot.ptr() ) {
Rotation3D rot3D(RotationZYX(rot.z(0),rot.y(0),rot.x(0)));
Position pos3D(pos.x(0),pos.y(0),pos.z(0));
tr = Transform3D(rot3D, pos3D);
}
else if ( pos.ptr() ) {
tr = Transform3D(Rotation3D(),Position(pos.x(0),pos.y(0),pos.z(0)));
}
else if ( rot.ptr() ) {
Rotation3D rot3D(RotationZYX(rot.z(0),rot.y(0),rot.x(0)));
tr = Transform3D(rot3D,Position());
}
Material support_mat = description.material(x_support.materialStr());
Volume support_vol(support_name, support_solid, support_mat);
support_vol.setVisAttributes(description.visAttributes(support_vis));
pv = assembly.placeVolume(support_vol, tr);
// pv = assembly.placeVolume(support_vol, Position(0, 0, support_zstart + support_length / 2));
}
for (xml_coll_t mi(x_det, _U(module)); mi; ++mi, ++m_id) { for (xml_coll_t mi(x_det, _U(module)); mi; ++mi, ++m_id) {
xml_comp_t x_mod = mi; xml_comp_t x_mod = mi;
string m_nam = x_mod.nameStr(); string m_nam = x_mod.nameStr();
...@@ -60,12 +122,15 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -60,12 +122,15 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
double x1 = trd.x1(); double x1 = trd.x1();
double x2 = trd.x2(); double x2 = trd.x2();
double z = trd.z(); double z = trd.z();
double y1, y2, total_thickness = 0.; double total_thickness = 0.;
xml_coll_t ci(x_mod, _U(module_component)); xml_coll_t ci(x_mod, _U(module_component));
for (ci.reset(), total_thickness = 0.0; ci; ++ci) for (ci.reset(), total_thickness = 0.0; ci; ++ci)
total_thickness += xml_comp_t(ci).thickness(); total_thickness += xml_comp_t(ci).thickness();
y1 = y2 = total_thickness / 2; double thickness_so_far = 0.0;
double thickness_sum = -total_thickness / 2.0;
double y1 = total_thickness / 2;
double y2 = total_thickness / 2;
Trapezoid m_solid(x1, x2, y1, y2, z); Trapezoid m_solid(x1, x2, y1, y2, z);
Volume m_volume(m_nam, m_solid, vacuum); Volume m_volume(m_nam, m_solid, vacuum);
m_volume.setVisAttributes(description.visAttributes(x_mod.visStr())); m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
...@@ -114,13 +179,38 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -114,13 +179,38 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
c_vol.setVisAttributes(description.visAttributes(c.visStr())); c_vol.setVisAttributes(description.visAttributes(c.visStr()));
pv = m_volume.placeVolume(c_vol, Position(0, posY + c_thick / 2, 0)); pv = m_volume.placeVolume(c_vol, Position(0, posY + c_thick / 2, 0));
if (c.isSensitive()) { if (c.isSensitive()) {
module_thicknesses[m_nam] = {thickness_so_far + c_thick/2.0, total_thickness-thickness_so_far - c_thick/2.0};
//std::cout << " adding sensitive volume" << c_name << "\n";
sdet.check(n_sensor > 2, "SiTrackerEndcap2::fromCompact: " + c_name + " Max of 2 modules allowed!"); sdet.check(n_sensor > 2, "SiTrackerEndcap2::fromCompact: " + c_name + " Max of 2 modules allowed!");
pv.addPhysVolID("sensor", n_sensor); pv.addPhysVolID("sensor", n_sensor);
sens.setType("tracker");
c_vol.setSensitiveDetector(sens); c_vol.setSensitiveDetector(sens);
sensitives[m_nam].push_back(pv); sensitives[m_nam].push_back(pv);
++n_sensor; ++n_sensor;
// -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
Vector3D u(0., 0., -1.);
Vector3D v(-1., 0., 0.);
Vector3D n(0., 1., 0.);
// Vector3D o( 0. , 0. , 0. ) ;
// compute the inner and outer thicknesses that need to be assigned to the tracking surface
// depending on wether the support is above or below the sensor
double inner_thickness = module_thicknesses[m_nam][0];
double outer_thickness = module_thicknesses[m_nam][1];
SurfaceType type(SurfaceType::Sensitive);
// if( isStripDetector )
// type.setProperty( SurfaceType::Measurement1D , true ) ;
VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
volplane_surfaces[m_nam].push_back(surf);
//--------------------------------------------
} }
posY += c_thick; posY += c_thick;
thickness_sum += c_thick;
thickness_so_far += c_thick;
} }
modules[m_nam] = m_volume; modules[m_nam] = m_volume;
} }
...@@ -152,18 +242,25 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -152,18 +242,25 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
if (reflect) { if (reflect) {
layer_pv = layer_pv =
assembly.placeVolume(layer_vol, Transform3D(RotationZYX(0.0, -M_PI, 0.0), Position(0, 0, -layer_center_z))); assembly.placeVolume(layer_vol, Transform3D(RotationZYX(0.0, -M_PI, 0.0), Position(0, 0, -layer_center_z)));
layer_pv.addPhysVolID("barrel", 3).addPhysVolID("layer", l_id); layer_pv.addPhysVolID("layer", l_id);
layer_name += "_N"; layer_name += "_N";
} else { } else {
layer_pv = assembly.placeVolume(layer_vol, Position(0, 0, layer_center_z)); layer_pv = assembly.placeVolume(layer_vol, Position(0, 0, layer_center_z));
layer_pv.addPhysVolID("barrel", 2).addPhysVolID("layer", l_id); layer_pv.addPhysVolID("layer", l_id);
layer_name += "_P"; layer_name += "_P";
} }
DetElement layer_element(sdet, layer_name, l_id); DetElement layer_element(sdet, layer_name, l_id);
layer_element.setPlacement(layer_pv); layer_element.setPlacement(layer_pv);
Acts::ActsExtension* layerExtension = new Acts::ActsExtension(); Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
layerExtension->addType("layer", "layer"); layerExtension->addType("sensitive disk", "layer");
//layerExtension->addType("axes", "definitions", "XZY"); //layerExtension->addType("axes", "definitions", "XZY");
//layerExtension->addType("sensitive disk", "layer");
//layerExtension->addType("axes", "definitions", "XZY");
for (xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
xml_comp_t x_layer_material = lmat;
xmlToProtoSurfaceMaterial(x_layer_material, *layerExtension, "layer_material");
}
layer_element.addExtension<Acts::ActsExtension>(layerExtension); layer_element.addExtension<Acts::ActsExtension>(layerExtension);
for (xml_coll_t ri(x_layer, _U(ring)); ri; ++ri) { for (xml_coll_t ri(x_layer, _U(ring)); ri; ++ri) {
...@@ -185,30 +282,34 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -185,30 +282,34 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
double y = -r * std::sin(phi); double y = -r * std::sin(phi);
if (!reflect) { if (!reflect) {
DetElement module(sdet, m_base + "_pos", det_id); DetElement module(layer_element, m_base + "_pos", det_id);
pv = layer_vol.placeVolume( pv = layer_vol.placeVolume(
m_vol, Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, zstart + dz))); m_vol, Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, zstart + dz)));
pv.addPhysVolID("barrel", 1).addPhysVolID("layer", l_id).addPhysVolID("module", mod_num); pv.addPhysVolID("module", mod_num);
module.setPlacement(pv); module.setPlacement(pv);
for (size_t ic = 0; ic < sensVols.size(); ++ic) { for (size_t ic = 0; ic < sensVols.size(); ++ic) {
PlacedVolume sens_pv = sensVols[ic]; PlacedVolume sens_pv = sensVols[ic];
DetElement comp_elt(module, sens_pv.volume().name(), mod_num); DetElement comp_elt(module, sens_pv.volume().name(), mod_num);
comp_elt.setPlacement(sens_pv); comp_elt.setPlacement(sens_pv);
Acts::ActsExtension* moduleExtension = new Acts::ActsExtension(); //std::cout << " adding ACTS extension" << "\n";
Acts::ActsExtension* moduleExtension = new Acts::ActsExtension("XZY");
comp_elt.addExtension<Acts::ActsExtension>(moduleExtension); comp_elt.addExtension<Acts::ActsExtension>(moduleExtension);
volSurfaceList(comp_elt)->push_back(volplane_surfaces[m_nam][ic]);
} }
} else { } else {
pv = layer_vol.placeVolume( pv = layer_vol.placeVolume(
m_vol, Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, -zstart - dz))); m_vol, Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, -zstart - dz)));
pv.addPhysVolID("barrel", 2).addPhysVolID("layer", l_id).addPhysVolID("module", mod_num); pv.addPhysVolID("module", mod_num);
DetElement r_module(sdet, m_base + "_neg", det_id); DetElement r_module(layer_element, m_base + "_neg", det_id);
r_module.setPlacement(pv); r_module.setPlacement(pv);
for (size_t ic = 0; ic < sensVols.size(); ++ic) { for (size_t ic = 0; ic < sensVols.size(); ++ic) {
PlacedVolume sens_pv = sensVols[ic]; PlacedVolume sens_pv = sensVols[ic];
DetElement comp_elt(r_module, sens_pv.volume().name(), mod_num); DetElement comp_elt(r_module, sens_pv.volume().name(), mod_num);
comp_elt.setPlacement(sens_pv); comp_elt.setPlacement(sens_pv);
Acts::ActsExtension* moduleExtension = new Acts::ActsExtension(); //std::cout << " adding ACTS extension" << "\n";
Acts::ActsExtension* moduleExtension = new Acts::ActsExtension("XZY");
comp_elt.addExtension<Acts::ActsExtension>(moduleExtension); comp_elt.addExtension<Acts::ActsExtension>(moduleExtension);
volSurfaceList(comp_elt)->push_back(volplane_surfaces[m_nam][ic]);
} }
} }
dz = -dz; dz = -dz;
...@@ -223,6 +324,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -223,6 +324,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
return sdet; return sdet;
} }
//@}
// clang-format off // clang-format off
DECLARE_DETELEMENT(refdet_TrapEndcapTracker, create_detector) DECLARE_DETELEMENT(athena_TrapEndcapTracker, create_detector)
DECLARE_DETELEMENT(refdet_GEMTrackerEndcap, create_detector) DECLARE_DETELEMENT(athena_GEMTrackerEndcap, create_detector)
src/allsilicon_support.cpp deleted 100644 → 0
View file @ 2e7fa037
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/Printout.h"
#include "DDRec/DetectorData.h"
#include "DDRec/Surface.h"
#include <XML/Helper.h>
//////////////////////////////////
// Support structure for ALl-silicon
//////////////////////////////////
using namespace std;
using namespace dd4hep;
// Info from
// https://github.com/reynier0611/g4lblvtx/blob/master/source/AllSi_vtx_serv_2lyr_Detector.cc
// TODO: this is quite incomplete, should probably wait for official word
// from he tracking WG
static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
{
xml_det_t x_det = e;
string detName = x_det.nameStr();
int detID = x_det.id();
bool reflect = x_det.reflect();
const int sign = reflect ? -1 : 1;
// second vertexing layer
std::vector<double> z_det = {15 * cm, 20 * cm};
std::vector<double> rin_l2 = {5.48 * cm, 14.8 * cm};
std::vector<double> rout_l2 = {0, 0};
// first vertexing layer
std::vector<double> rin_l1 = {3.30 * cm, 14.36 * cm};
std::vector<double> rout_l1 = {0, 0};
const int nzplanes_l2 = z_det.size();
const int nzplanes_l1 = z_det.size();
for (int i = 0; i < nzplanes_l2; i++) {
rout_l2[i] = rin_l2[i] + 0.44;
z_det[i] *= sign / abs(sign);
}
for (int i = 0; i < nzplanes_l1; i++) {
rout_l1[i] = rin_l1[i] + 0.44;
z_det[i] *= sign / abs(sign);
}
// mother volume
std::vector<double> rin_mo = rin_l1;
std::vector<double> rout_mo = rout_l2;
DetElement det(detName, detID);
Material Vacuum = desc.material("Vacuum");
Polycone empty_cone("empty_cone", 0.0, 360 * degree, z_det, rin_mo, rout_mo);
Volume detVol("empty_cone", empty_cone, Vacuum);
detVol.setVisAttributes(desc.invisible());
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(0.0, 0.0, 0.0), Position(0.0, 0.0, 0.));
PlacedVolume detPV = motherVol.placeVolume(detVol, tr);
detPV.addPhysVolID("system", detID);
detPV.addPhysVolID("barrel", 1);
det.setPlacement(detPV);
Material Al = desc.material("Al");
Material Graphite = desc.material("Graphite");
// cb_DIRC_bars_Logic.setVisAttributes(desc.visAttributes(x_det.visStr()));
Polycone polycone_l2("polycone_l2", 0, 360 * degree, z_det, rin_l2, rout_l2);
Volume logical_l2("polycone_l2_logic", polycone_l2, Al);
logical_l2.setVisAttributes(desc.visAttributes(x_det.visStr()));
detVol.placeVolume(logical_l2, tr);
Polycone polycone_l1("polycone_l1", 0, 360 * degree, z_det, rin_l1, rout_l1);
Volume logical_l1("polycone_l1_logic", polycone_l1, Al);
logical_l1.setVisAttributes(desc.visAttributes(x_det.visStr()));
detVol.placeVolume(logical_l1, tr);
return det;
}
DECLARE_DETELEMENT(allsilicon_support, createDetector)
src/cb_CTD_Si.cpp deleted 100644 → 0
View file @ 2e7fa037
#include <XML/Helper.h>
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
//////////////////////////////////
// Central Barrel Tracker Silicon
//////////////////////////////////
using namespace std;
using namespace dd4hep;
static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
{
xml_det_t x_det = e;
string detName = x_det.nameStr();
int detID = x_det.id();
xml_dim_t dim = x_det.dimensions();
double RIn = dim.rmin();
double ROut = dim.rmax();
double SizeZ = dim.length();
double SizeZCut = dim.zmax();
double SiLayerGap = dim.gap();
Material Vacuum = desc.material("Vacuum");
// Create Global Volume
Tube cb_CTD_GVol_Solid(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
Volume detVol("cb_CTD_GVol_Logic", cb_CTD_GVol_Solid, Vacuum);
detVol.setVisAttributes(desc.visAttributes(x_det.visStr()));
// Construct Silicon Layers
xml_comp_t x_layer = x_det.child(_U(layer));
const int repeat = x_layer.repeat();
xml_comp_t x_slice = x_layer.child(_U(slice));
Material slice_mat = desc.material(x_slice.materialStr());
double layerRIn[100];
double layerROut[100];
// Loop over layers
for(int i = 0; i < repeat; i++) {
layerRIn[i] = RIn + (SiLayerGap * i);
layerROut[i] = RIn + (0.01 + SiLayerGap * i);
if (layerROut[i] > ROut)
continue;
string logic_layer_name = detName + _toString(i, "_Logic_lay_%d");
if (i==7){logic_layer_name = detName + _toString(20, "_Logic_lay_%d");}
Volume layerVol(logic_layer_name,Tube(layerRIn[i], layerROut[i], SizeZ / 2.0, 0.0, 360.0 * deg), slice_mat);
layerVol.setVisAttributes(desc,x_layer.visStr());
sens.setType("tracker");
layerVol.setSensitiveDetector(sens);
Position layer_pos = Position(0.0, 0.0, 0.0);
PlacedVolume layerPV = detVol.placeVolume(layerVol, layer_pos);
layerPV.addPhysVolID("layer", i+1);
}
DetElement det(detName, detID);
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(0.0, 0.0, 0.0), Position(0.0, 0.0, 0.0));
PlacedVolume detPV = motherVol.placeVolume(detVol, tr);
detPV.addPhysVolID("system", detID);
detPV.addPhysVolID("barrel", 1);
det.setPlacement(detPV);
return det;
}
DECLARE_DETELEMENT(cb_CTD_Si, createDetector)
src/cb_DIRC.cpp deleted 100644 → 0
View file @ 2e7fa037
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/Printout.h"
#include "DDRec/DetectorData.h"
#include "DDRec/Surface.h"
#include <XML/Helper.h>
//////////////////////////////////
// Central Barrel DIRC
//////////////////////////////////
using namespace std;
using namespace dd4hep;
static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
{
xml_det_t x_det = e;
string detName = x_det.nameStr();
int detID = x_det.id();
xml_dim_t dim = x_det.dimensions();
xml_dim_t pos = x_det.position();
double RIn = dim.rmin();
double ROut = dim.rmax();
double SizeZ = dim.length();
Material Vacuum = desc.material("Vacuum");
Tube cb_DIRC_Barrel_GVol_Solid(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
Volume detVol("cb_DIRC_GVol_Solid_Logic", cb_DIRC_Barrel_GVol_Solid, Vacuum);
detVol.setVisAttributes(desc.invisible());
DetElement det(detName, detID);
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(0.0, 0.0, 0.0), Position(0.0, 0.0, pos.z()));
PlacedVolume detPV = motherVol.placeVolume(detVol, tr);
detPV.addPhysVolID("system", detID);
detPV.addPhysVolID("barrel", 1);
det.setPlacement(detPV);
//////////////////
// DIRC Bars
//////////////////
double dR = 83.65 * cm;
double cb_DIRC_bars_DZ = SizeZ;
double cb_DIRC_bars_DY = 42. * cm;
double cb_DIRC_bars_DX = 1.7 * cm;
double myL = 2 * M_PI * dR;
int NUM = myL / cb_DIRC_bars_DY;
double cb_DIRC_bars_deltaphi = 2 * 3.1415926 / NUM;
Material cb_DIRC_bars_Material = desc.material("Quartz");
Box cb_DIRC_bars_Solid(cb_DIRC_bars_DX / 2., cb_DIRC_bars_DY / 2., cb_DIRC_bars_DZ / 2.);
Volume cb_DIRC_bars_Logic("cb_DIRC_bars_Logix", cb_DIRC_bars_Solid, cb_DIRC_bars_Material);
cb_DIRC_bars_Logic.setVisAttributes(desc.visAttributes(x_det.visStr()));
sens.setType("photoncounter");
cb_DIRC_bars_Logic.setSensitiveDetector(sens);
for (int ia = 0; ia < NUM; ia++) {
double phi = (ia * (cb_DIRC_bars_deltaphi));
double x = -dR * cos(phi);
double y = -dR * sin(phi);
Transform3D tr(RotationZ(cb_DIRC_bars_deltaphi * ia), Position(x, y, 0));
PlacedVolume barPV = detVol.placeVolume(cb_DIRC_bars_Logic, tr);
barPV.addPhysVolID("module", ia);
}
return det;
}
DECLARE_DETELEMENT(cb_DIRC, createDetector)
src/cb_VTX_Barrel.cpp deleted 100644 → 0
View file @ 2e7fa037
#include <XML/Helper.h>
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
//////////////////////////////////
// Central Barrel Vertex Detector
//////////////////////////////////
using namespace std;
using namespace dd4hep;
static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
{
xml_det_t x_det = e;
string detName = x_det.nameStr();
int detID = x_det.id();
xml_dim_t dim = x_det.dimensions();
double RIn = dim.rmin();
double ROut = dim.rmax();
double SizeZ = dim.length();
xml_dim_t pos = x_det.position();
Material Vacuum = desc.material("Vacuum");
// Create Global Volume
Tube cb_VTX_Barrel_GVol_Solid(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
Volume detVol("cb_VTX_Barrel_GVol_Logic", cb_VTX_Barrel_GVol_Solid, Vacuum);
detVol.setVisAttributes(desc.visAttributes(x_det.visStr()));
//////////////////
// Barrel Ladder
//////////////////
xml_comp_t x_layer = x_det.child(_U(layer));
const int repeat = x_layer.repeat();
xml_comp_t x_slice = x_layer.child(_U(slice));
Material slice_mat = desc.material(x_slice.materialStr());
double x = 0.0 * cm;
double y = 0.0 * cm;
double z = 0.0 * cm;
int FDIV = 0;
double dR = 0.0;
double length = 0.0;
double phi = 0.0;
// Ladder Layer Parameters
double lay_Dx[6];
double lay_Dy[6];
double lay_Dz[6];
double lay_Rin[6];
lay_Dx[0] = 0.050 * mm; lay_Dy[0] = 1.0 * cm; lay_Dz[0] = 10.0 * cm; lay_Rin[0] = 3.5 * cm;
lay_Dx[1] = 0.050 * mm; lay_Dy[1] = 1.0 * cm; lay_Dz[1] = 11.0 * cm; lay_Rin[1] = 4.5 * cm;
lay_Dx[2] = 0.150 * mm; lay_Dy[2] = 2.0 * cm; lay_Dz[2] = 18.0 * cm; lay_Rin[2] = 6.5 * cm;
lay_Dx[3] = 0.150 * mm; lay_Dy[3] = 2.0 * cm; lay_Dz[3] = 24.0 * cm; lay_Rin[3] = 10.5 * cm;
lay_Dx[4] = 0.150 * mm; lay_Dy[4] = 3.0 * cm; lay_Dz[4] = 36.0 * cm; lay_Rin[4] = 13.5 * cm;
lay_Dx[5] = 0.150 * mm; lay_Dy[5] = 3.0 * cm; lay_Dz[5] = 48.0 * cm; lay_Rin[5] = 15.5 * cm;
int i_layer = 0;
int i_module = 0;
// Loop over layers
for(int i = 0; i < repeat; i++) {
double cb_VTX_Barrel_ladder_DZ = lay_Dz[i];
double cb_VTX_Barrel_ladder_DY = lay_Dy[i];
double cb_VTX_Barrel_ladder_Thickness = lay_Dx[i];
dR = lay_Rin[i];
length = 2.0 * 3.1415 * dR;
int laddersCount = length / cb_VTX_Barrel_ladder_DY;
for (int i = 0; i < 2; i++) {
double LN = cb_VTX_Barrel_ladder_DY * laddersCount;
double LN1 = cb_VTX_Barrel_ladder_DY * (laddersCount + 1.0 + i);
if (LN/LN1 > 0.8)
laddersCount = laddersCount + 1;
}
double cb_VTX_Barrel_ladder_deltaphi = 2.0 * 3.1415926 / laddersCount;
string ladderBoxName = detName + _toString(i, "_ladder_Solid_%d");
string ladderName = detName + _toString(i, "_ladder_Logic_%d");
Volume ladderVol(ladderName, Box(cb_VTX_Barrel_ladder_Thickness * 0.5, cb_VTX_Barrel_ladder_DY * 0.5, cb_VTX_Barrel_ladder_DZ * 0.5), slice_mat);
ladderVol.setVisAttributes(desc,x_layer.visStr());
sens.setType("tracker");
ladderVol.setSensitiveDetector(sens);
i_layer++;
for (int ia = 0; ia < laddersCount; ia++) {
phi = (ia * (cb_VTX_Barrel_ladder_deltaphi));
x = - dR * cos(phi);
y = - dR * sin(phi);
RotationZYX ladder_rot = RotationZYX(cb_VTX_Barrel_ladder_deltaphi * ia, 0.0, 0.0);
Position ladder_pos = Position(x, y, z);
string ladderName = detName + _toString(i, "_ladder_Phys_%d") + _toString(ia, "_%d");
PlacedVolume ladderPV = detVol.placeVolume(ladderVol, Transform3D(ladder_rot, ladder_pos));
i_module++;
ladderPV.addPhysVolID("layer", i_layer).addPhysVolID("module", i_module);
}
}
// TODO: Pixels
DetElement det(detName, detID);
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(0.0, 0.0, 0.0), Position(0.0, 0.0, pos.z()));
PlacedVolume detPV = motherVol.placeVolume(detVol, tr);
detPV.addPhysVolID("system", detID);
detPV.addPhysVolID("barrel", 1);
det.setPlacement(detPV);
return det;
}
DECLARE_DETELEMENT(cb_VTX_Barrel, createDetector)
src/ce_MRICH.cpp deleted 100644 → 0
View file @ 2e7fa037
#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"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, SensitiveDetector &sens);
// 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();
// xml::Component rads = detElem.child(_Unicode(radiator));
auto rmin = dims.rmin();
auto rmax = dims.rmax();
auto length = dims.length();
auto z0 = dims.z();
auto gasMat = desc.material("AirOptical");
// detector envelope
Tube envShape(rmin, rmax, length/2., 0., 2*M_PI);
Volume envVol("ce_MRICH_GVol", envShape, gasMat);
envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
// modules
addModules(envVol, detElem, desc, sens);
// place envelope
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0));
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, 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 = desc.material(mods.materialStr());
auto gasMat = desc.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(desc.visAttributes(mods.visStr()));
sens.setType("photoncounter");
modVol.setSensitiveDetector(sens);
// place modules in the sectors (disk)
auto points = athena::geo::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(refdet_ce_MRICH, createDetector)
src/BarrelCalorimeterHybrid_geo.cpp → src/deprecated/BarrelCalorimeterHybrid_geo.cpp
View file @ 1008a184
...@@ -29,37 +29,41 @@ using namespace dd4hep::detail; ...@@ -29,37 +29,41 @@ using namespace dd4hep::detail;
typedef ROOT::Math::XYPoint Point; typedef ROOT::Math::XYPoint Point;
// Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system // Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
vector<Point> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi, double spacing_tol = 1e-2) { vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi, double spacing_tol = 1e-2) {
// z_spacing - distance between fiber layers in z // z_spacing - distance between fiber layers in z
// x_spacing - distance between fiber centers in x // x_spacing - distance between fiber centers in x
// x - half-length of the shorter (bottom) base of the trapezoid // x - half-length of the shorter (bottom) base of the trapezoid
// z - height of the trapezoid // z - height of the trapezoid
// phi - angle between z and trapezoid arm // phi - angle between z and trapezoid arm
vector<Point> positions; vector<vector<Point>> positions;
int z_layers = floor((z/2-radius-spacing_tol)/z_spacing); // number of layers that fit in z/2 int z_layers = floor((z/2-radius-spacing_tol)/z_spacing); // number of layers that fit in z/2
double z_pos = 0.; double z_pos = 0.;
double x_pos = 0.; double x_pos = 0.;
for(int l = -z_layers; l < z_layers+1; l++) { for(int l = -z_layers; l < z_layers+1; l++) {
vector<Point> xline;
z_pos = l*z_spacing;
double x_max = x + (z/2. + z_pos)*tan(phi) - spacing_tol; // calculate max x at particular z_pos
(l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing/2; // account for spacing/2 shift
z_pos = l*z_spacing; while(x_pos < (x_max - radius)) {
double x_max = x + (z/2. + z_pos)*tan(phi) - spacing_tol; // calculate max x at particular z_pos xline.push_back(Point(x_pos, z_pos));
(l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing/2; // account for spacing/2 shift if(x_pos != 0.) xline.push_back(Point(-x_pos, z_pos)); // using symmetry around x=0
x_pos += x_spacing;
while(x_pos < (x_max - radius)) {
positions.push_back(Point(x_pos,z_pos));
if(x_pos != 0.) positions.push_back(Point(-x_pos,z_pos)); // using symmetry around x=0
x_pos += x_spacing;
}
} }
// Sort fiber IDs for a better organization
return positions; sort(xline.begin(), xline.end(), [](const Point &p1, const Point &p2) {
return p1.x() < p2.x();
});
positions.emplace_back(std::move(xline));
}
return positions;
} }
// Calculate number of divisions for the readout grid for the fiber layers // Calculate number of divisions for the readout grid for the fiber layers
std::pair<int, int> getNdivisions(double x, double z, double dx, double dz){ std::pair<int, int> getNdivisions(double x, double z, double dx, double dz) {
// x and z defined as in vector<Point> fiberPositions // x and z defined as in vector<Point> fiberPositions
// dx, dz - size of the grid in x and z we want to get close to with the polygons // dx, dz - size of the grid in x and z we want to get close to with the polygons
// See also descripltion when the function is called // See also descripltion when the function is called
...@@ -302,14 +306,6 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -302,14 +306,6 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
std::string f_id_fiber = getAttrOrDefault(x_fiber, _Unicode(identifier_fiber), "fiber"); std::string f_id_fiber = getAttrOrDefault(x_fiber, _Unicode(identifier_fiber), "fiber");
// Calculate fiber positions inside the slice // Calculate fiber positions inside the slice
vector<Point> f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick-tolerance, hphi);
// Sort fiber IDs fo better organization
sort(f_pos.begin(), f_pos.end(),
[](const Point &p1, const Point &p2) {
if (p1.y() == p2.y()) { return p1.x() < p2.x(); }
return p1.y() < p2.y();
});
Tube f_tube(0, f_radius, stave_z-tolerance); Tube f_tube(0, f_radius, stave_z-tolerance);
// Set up the readout grid for the fiber layers // Set up the readout grid for the fiber layers
...@@ -325,42 +321,44 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -325,42 +321,44 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
vector<tuple<int, Point, Point, Point, Point>> grid_vtx = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick-tolerance, hphi); vector<tuple<int, Point, Point, Point, Point>> grid_vtx = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick-tolerance, hphi);
vector<int> f_id_count(grid_div.first*grid_div.second,0); vector<int> f_id_count(grid_div.first*grid_div.second,0);
for (auto &p : f_pos) { auto f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick-tolerance, hphi);
int f_grid_id = -1; for (auto &line : f_pos) {
int f_id = -1; for (auto &p : line) {
// Check to which grid fiber belongs to int f_grid_id = -1;
for (auto &poly_vtx : grid_vtx) { int f_id = -1;
auto [grid_id, vtx_a, vtx_b, vtx_c, vtx_d] = poly_vtx; // Check to which grid fiber belongs to
double poly_x[4] = {vtx_a.x(), vtx_b.x(), vtx_c.x(), vtx_d.x()}; for (auto &poly_vtx : grid_vtx) {
double poly_y[4] = {vtx_a.y(), vtx_b.y(), vtx_c.y(), vtx_d.y()}; auto [grid_id, vtx_a, vtx_b, vtx_c, vtx_d] = poly_vtx;
double f_xy[2] = {p.x(), p.y()}; double poly_x[4] = {vtx_a.x(), vtx_b.x(), vtx_c.x(), vtx_d.x()};
double poly_y[4] = {vtx_a.y(), vtx_b.y(), vtx_c.y(), vtx_d.y()};
TGeoPolygon poly(4); double f_xy[2] = {p.x(), p.y()};
poly.SetXY(poly_x,poly_y);
poly.FinishPolygon(); TGeoPolygon poly(4);
poly.SetXY(poly_x,poly_y);
if(poly.Contains(f_xy)) { poly.FinishPolygon();
f_grid_id = grid_id;
f_id = f_id_count[grid_id]; if(poly.Contains(f_xy)) {
f_id_count[grid_id]++; f_grid_id = grid_id;
f_id = f_id_count[grid_id];
f_id_count[grid_id]++;
}
} }
}
string f_name = "fiber" + to_string(f_grid_id) + "_" + to_string(f_id);
Volume f_vol(f_name, f_tube, description.material(x_fiber.materialStr()));
DetElement fiber(slice, f_name, det_id);
if ( x_fiber.isSensitive() ) {
f_vol.setSensitiveDetector(sens);
}
fiber.setAttributes(description,f_vol,x_fiber.regionStr(),x_fiber.limitsStr(),x_fiber.visStr());
// Fiber placement string f_name = "fiber" + to_string(f_grid_id) + "_" + to_string(f_id);
Transform3D f_tr(RotationZYX(0,0,M_PI*0.5),Position(p.x(), 0 ,p.y())); Volume f_vol(f_name, f_tube, description.material(x_fiber.materialStr()));
PlacedVolume fiber_phv = s_vol.placeVolume(f_vol, f_tr); DetElement fiber(slice, f_name, det_id);
fiber_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1); if ( x_fiber.isSensitive() ) {
fiber.setPlacement(fiber_phv); f_vol.setSensitiveDetector(sens);
}
fiber.setAttributes(description,f_vol,x_fiber.regionStr(),x_fiber.limitsStr(),x_fiber.visStr());
} // Fiber placement
Transform3D f_tr(RotationZYX(0,0,M_PI*0.5),Position(p.x(), 0 ,p.y()));
PlacedVolume fiber_phv = s_vol.placeVolume(f_vol, f_tr);
fiber_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1);
fiber.setPlacement(fiber_phv);
}
}
} }
if ( x_slice.isSensitive() ) { if ( x_slice.isSensitive() ) {
......
views/generate_prim_files.yml 0 → 100644
View file @ 1008a184
view_prim:detector_only:
extends: .views
stage: test
script:
- ./bin/generate_prim_file -o ${LOCAL_DATA_PATH} -D -t detector_view
- ls -lrth && ls -lrth ${LOCAL_DATA_PATH}
view_prim:ev001:
extends: .views
stage: test
rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"'
script:
- ./bin/generate_prim_file -o ${LOCAL_DATA_PATH} -t view_ev001 -s 1
view_prim:ev002:
extends: .views
stage: test
rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"'
script:
- ./bin/generate_prim_file -o ${LOCAL_DATA_PATH} -t view_ev002 -s 2
view_prim:ev003:
extends: .views
stage: test
rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"'
script:
- ./bin/generate_prim_file -o ${LOCAL_DATA_PATH} -t view_ev003 -s 3
view_prim:ev004:
extends: .views
stage: test
rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"'
script:
- ./bin/generate_prim_file -o ${LOCAL_DATA_PATH} -t view_ev004 -s 4
view_prim:calorimeters:
extends: .views
stage: test
script:
- cp "compact/subsystem_views/calorimeters.xml" "${DETECTOR_PATH}/."
- ./bin/generate_prim_file -c ${DETECTOR_PATH}/calorimeters.xml -o ${LOCAL_DATA_PATH} -D -t calorimeters_view
- ls -lrth && ls -lrth ${LOCAL_DATA_PATH}
view_prim:calorimeters_ev001:
extends: .views
stage: test
rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"'
script:
- cp "compact/subsystem_views/calorimeters.xml" "${DETECTOR_PATH}/."
- ./bin/generate_prim_file -c ${DETECTOR_PATH}/calorimeters.xml -o ${LOCAL_DATA_PATH} -t calorimeters_view_ev001 -s 1
view_prim:calorimeters_ev002:
extends: .views
stage: test
rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"'
script:
- cp "compact/subsystem_views/calorimeters.xml" "${DETECTOR_PATH}/."
- ./bin/generate_prim_file -c ${DETECTOR_PATH}/calorimeters.xml -o ${LOCAL_DATA_PATH} -t calorimeters_view_ev002 -s 2
views/view1.yml
View file @ 1008a184
dawn_view_01:detector: dawn_view_01:detector:
extends: .views extends: .views
needs:
- job: view_prim:detector_only
optional: false
script: script:
- ./bin/make_dawn_views -t view01 -d scripts/view1 -D - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/detector_view.prim -t view01 -d scripts/view1 -D
dawn_view_01:ev001: dawn_view_01:ev001:
extends: .views extends: .views
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
needs:
- job: view_prim:ev001
optional: true
script: script:
- ./bin/make_dawn_views -t view01_ev001 -d scripts/view1 -s 1 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/view_ev001.prim -t view01_ev001 -d scripts/view1 -s 1
dawn_view_01:ev002: dawn_view_01:ev002:
extends: .views extends: .views
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
needs:
- job: view_prim:ev002
optional: true
script: script:
- ./bin/make_dawn_views -t view01_ev002 -d scripts/view1 -s 2 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/view_ev002.prim -t view01_ev002 -d scripts/view1 -s 2
view_01: view_01:
stage: test stage: collect
rules:
- if: '$CI_SERVER_HOST == "eicweb.phy.anl.gov"'
needs: needs:
- job: dawn_view_01:detector - job: dawn_view_01:detector
optional: false optional: false
......
views/view11.yml
View file @ 1008a184
dawn_view_11:detector: dawn_view_11:detector:
extends: .views extends: .views
needs:
- job: view_prim:detector_only
optional: false
script: script:
- ./bin/make_dawn_views -t view11 -d scripts/view11 -D - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/detector_view.prim -t view11 -d scripts/view11 -D
dawn_view_11:ev000: dawn_view_11:ev000:
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
extends: .views extends: .views
needs:
- job: view_prim:ev001
optional: true
script: script:
- ./bin/make_dawn_views -t view11 -d scripts/view11 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/view_ev001.prim -t view11 -d scripts/view11
dawn_view_11:ev001: dawn_view_11:ev001:
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
extends: .views extends: .views
needs:
- job: view_prim:ev001
optional: true
script: script:
- ./bin/make_dawn_views -t view11 -d scripts/view11 -s 1 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/view_ev002.prim -t view11 -d scripts/view11 -s 1
dawn_view_11:ev002: dawn_view_11:ev002:
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
extends: .views extends: .views
needs:
- job: view_prim:ev002
optional: true
script: script:
- ./bin/make_dawn_views -t view11 -d scripts/view11 -s 2 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH}/view_ev003.prim -t view11 -d scripts/view11 -s 2
dawn_view_11:ev003: dawn_view_11:ev003:
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
extends: .views extends: .views
needs:
- job: view_prim:ev003
optional: true
script: script:
- ./bin/make_dawn_views -t view11 -d scripts/view11 -s 3 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH} -t view11 -d scripts/view11 -s 3
dawn_view_11:ev004: dawn_view_11:ev004:
rules: rules:
- if: '$DETECTOR_EVENT_VIEWS == "ON"' - if: '$DETECTOR_EVENT_VIEWS == "ON"'
extends: .views extends: .views
needs:
- job: view_prim:ev004
optional: true
script: script:
- ./bin/make_dawn_views -t view11 -d scripts/view11 -s 4 - ./bin/make_dawn_views -i ${LOCAL_DATA_PATH} -t view11 -d scripts/view11 -s 4
view_11: view_11:
stage: test stage: collect
rules: rules:
- if: '$CI_SERVER_HOST == "eicweb.phy.anl.gov"' - if: '$CI_SERVER_HOST == "eicweb.phy.anl.gov"'
needs: needs:
- job: compile
optional: false
- job: dawn_view_11:detector - job: dawn_view_11:detector
optional: false optional: false
- job: dawn_view_11:ev001 - job: dawn_view_11:ev001
......