g4simulation/g4micromegas/PHG4MicromegasDetector.cc

0001 /*!
0002  * \file PHG4MicromegasDetector.cc
0003  * \brief strongly inspired by code from Qinhua Huang <qinhua.huang@cea.fr>
0004  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0005  */
0006
0007 #include "PHG4MicromegasDetector.h"
0008
0009 #include "PHG4MicromegasDisplayAction.h"
0010 #include "PHG4MicromegasSurvey.h"
0011
0012 #include <phparameter/PHParameters.h>
0013
0014 #include <g4detectors/PHG4CylinderGeomContainer.h>
0015
0016 #include <g4main/PHG4Detector.h>
0017 #include <g4main/PHG4Subsystem.h>
0018
0019 #include <micromegas/CylinderGeomMicromegas.h>
0020
0021 #include <phool/PHCompositeNode.h>
0022 #include <phool/PHIODataNode.h>
0023 #include <phool/PHNode.h>          // for PHNode
0024 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0025 #include <phool/PHObject.h>        // for PHObject
0026 #include <phool/getClass.h>
0027 #include <phool/recoConsts.h>
0028
0029 #include <Geant4/G4Box.hh>
0030 #include <Geant4/G4Color.hh>
0031 #include <Geant4/G4LogicalVolume.hh>
0032 #include <Geant4/G4Material.hh>
0033 #include <Geant4/G4PVPlacement.hh>
0034 #include <Geant4/G4String.hh>  // for G4String
0035 #include <Geant4/G4SystemOfUnits.hh>
0036 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0037 #include <Geant4/G4Tubs.hh>
0038 #include <Geant4/G4Types.hh>            // for G4double
0039 #include <Geant4/G4VPhysicalVolume.hh>  // for G4VPhysicalVolume
0040 #include <Geant4/G4VSolid.hh>           // for G4VSolid
0041
0042 #include <cmath>
0043 #include <iostream>
0044 #include <numeric>
0045 #include <tuple>    // for make_tuple, tuple
0046 #include <utility>  // for pair, make_pair
0047 #include <vector>   // for vector
0048
0049 namespace
0050 {
0051   template <class T>
0052   constexpr T square(const T& x)
0053   {
0054     return x * x;
0055   }
0056   template <class T>
0057   inline T get_r(const T& x, const T& y)
0058   {
0059     return std::sqrt(square(x) + square(y));
0060   }
0061 }  // namespace
0062
0063 //____________________________________________________________________________..
0064 PHG4MicromegasDetector::PHG4MicromegasDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, PHParameters* parameters, const std::string& dnam)
0065   : PHG4Detector(subsys, Node, dnam)
0066   , m_DisplayAction(dynamic_cast<PHG4MicromegasDisplayAction*>(subsys->GetDisplayAction()))
0067   , m_Params(parameters)
0068   , m_ActiveFlag(m_Params->get_int_param("active"))
0069   , m_SupportActiveFlag(m_Params->get_int_param("supportactive"))
0070 {
0071   setup_tiles();
0072 }
0073
0074 //_______________________________________________________________
0075 int PHG4MicromegasDetector::IsInDetector(G4VPhysicalVolume* volume) const
0076 {
0077   if (m_ActiveFlag)
0078   {
0079     if (m_activeVolumes.contains(volume))
0080     {
0081       return 1;
0082     }
0083   }
0084   if (m_SupportActiveFlag)
0085   {
0086     if (m_passiveVolumes.contains(volume))
0087     {
0088       return -1;
0089     }
0090   }
0091   return 0;
0092 }
0093
0094 //_______________________________________________________________
0095 int PHG4MicromegasDetector::get_layer(G4VPhysicalVolume* volume) const
0096 {
0097   const auto iter = m_activeVolumes.find(volume);
0098   return iter == m_activeVolumes.end() ? -1 : iter->second;
0099 }
0100
0101 //_______________________________________________________________
0102 int PHG4MicromegasDetector::get_tileid(G4VPhysicalVolume* volume) const
0103 {
0104   const auto iter = m_tiles_map.find(volume);
0105   return iter == m_tiles_map.end() ? -1 : iter->second;
0106 }
0107
0108 //_______________________________________________________________
0109 void PHG4MicromegasDetector::ConstructMe(G4LogicalVolume* logicWorld)
0110 {
0111   create_materials();
0112   construct_micromegas(logicWorld);
0113   add_geometry_node();
0114 }
0115
0116 //_______________________________________________________________
0117 void PHG4MicromegasDetector::Print(const std::string& what) const
0118 {
0119   std::cout << "PHG4Micromegas Detector:" << std::endl;
0120   if (what == "ALL" || what == "VOLUME")
0121   {
0122     std::cout << "Version 0.1" << std::endl;
0123     std::cout << "Parameters:" << std::endl;
0124     m_Params->Print();
0125   }
0126   return;
0127 }
0128
0129 //_______________________________________________________________
0130 void PHG4MicromegasDetector::setup_tiles()
0131 {
0132   // TODO: replace with more realistic description from latest engineering drawings
0133   m_tiles.clear();
0134
0135   // tile dimensions
0136   /* they correspond to the total micromegas PCB size. They must match the definitions in construct_micromegas_tile */
0137   static constexpr double tile_length = 54.2;  // cm
0138   static constexpr double tile_width = 31.6;   // cm
0139
0140   {
0141     // bottom most sector 3pi/2 has 4 modules
0142     static constexpr double phi = 3. * M_PI / 2;
0143
0144     // tiles z position from CAD model (T-POT DETECTOR ASSEMBLY 7-13-2022-w-new-trays-modules)
0145     for (const double& tile_z : {-84.6, -28.2, 28.2, 84.6})
0146     {
0147       m_tiles.emplace_back(phi, tile_z, tile_width / CylinderGeomMicromegas::reference_radius, tile_length);
0148     }
0149   }
0150
0151   {
0152     // neighbor sectors have two modules, separated by 10cm
0153     for (const double& phi : {4. * M_PI / 3, 5. * M_PI / 3})
0154     {
0155       // tiles z position from CAD model (T-POT DETECTOR ASSEMBLY 7-13-2022-w-new-trays-modules)
0156       for (const double& tile_z : {-37.1, 37.1})
0157       {
0158         m_tiles.emplace_back(phi, tile_z, tile_width / CylinderGeomMicromegas::reference_radius, tile_length);
0159       }
0160     }
0161   }
0162 }
0163
0164 //_______________________________________________________________
0165 void PHG4MicromegasDetector::create_materials() const
0166 {
0167   // get the list of NIST materials
0168   // ---------------------------------
0169   auto *G4_N = GetDetectorMaterial("G4_N");
0170   auto *G4_O = GetDetectorMaterial("G4_O");
0171   auto *G4_C = GetDetectorMaterial("G4_C");
0172   auto *G4_H = GetDetectorMaterial("G4_H");
0173   auto *G4_Si = GetDetectorMaterial("G4_Si");
0174   auto *G4_Ar = GetDetectorMaterial("G4_Ar");
0175   auto *G4_Cr = GetDetectorMaterial("G4_Cr");
0176   auto *G4_Fe = GetDetectorMaterial("G4_Fe");
0177   auto *G4_Mn = GetDetectorMaterial("G4_Mn");
0178   auto *G4_Ni = GetDetectorMaterial("G4_Ni");
0179   auto *G4_Cu = GetDetectorMaterial("G4_Cu");
0180
0181   // combine elements
0182   // ----------------
0183   static const G4double temperature = 298.15 * kelvin;
0184   static const G4double pressure = 1. * atmosphere;
0185
0186   // FR4
0187   if (!GetDetectorMaterial("mmg_FR4", false))
0188   {
0189     auto *mmg_FR4 = new G4Material("mmg_FR4", 1.860 * g / cm3, 4, kStateSolid);
0190     mmg_FR4->AddMaterial(G4_C, 0.43550);
0191     mmg_FR4->AddMaterial(G4_H, 0.03650);
0192     mmg_FR4->AddMaterial(G4_O, 0.28120);
0193     mmg_FR4->AddMaterial(G4_Si, 0.24680);
0194   }
0195
0196   // Kapton
0197   if (!GetDetectorMaterial("mmg_Kapton", false))
0198   {
0199     auto *mmg_Kapton = new G4Material("mmg_Kapton", 1.420 * g / cm3, 4, kStateSolid);
0200     mmg_Kapton->AddMaterial(G4_C, 0.6911330);
0201     mmg_Kapton->AddMaterial(G4_H, 0.0263620);
0202     mmg_Kapton->AddMaterial(G4_N, 0.0732700);
0203     mmg_Kapton->AddMaterial(G4_O, 0.2092350);
0204   }
0205
0206   // MMgas
0207   if (!GetDetectorMaterial("mmg_Gas", false))
0208   {
0209     auto *mmg_Gas = new G4Material("mmg_Gas", 0.00170335 * g / cm3, 3, kStateGas, temperature, pressure);
0210     mmg_Gas->AddMaterial(G4_Ar, 0.900);
0211     mmg_Gas->AddMaterial(G4_C, 0.0826586);
0212     mmg_Gas->AddMaterial(G4_H, 0.0173414);
0213   }
0214
0215   // MMMesh
0216   if (!GetDetectorMaterial("mmg_Mesh", false))
0217   {
0218     auto *mmg_Mesh = new G4Material("mmg_Mesh", 2.8548 * g / cm3, 5, kStateSolid);
0219     mmg_Mesh->AddMaterial(G4_Cr, 0.1900);
0220     mmg_Mesh->AddMaterial(G4_Fe, 0.6800);
0221     mmg_Mesh->AddMaterial(G4_Mn, 0.0200);
0222     mmg_Mesh->AddMaterial(G4_Ni, 0.1000);
0223     mmg_Mesh->AddMaterial(G4_Si, 0.0100);
0224   }
0225
0226   // MMStrips
0227   if (!GetDetectorMaterial("mmg_Strips", false))
0228   {
0229     new G4Material("mmg_Strips", 5.248414 * g / cm3, G4_Cu, kStateSolid);
0230   }
0231
0232   // MMResistivePaste
0233   if (!GetDetectorMaterial("mmg_ResistPaste", false))
0234   {
0235     new G4Material("mmg_ResistPaste", 0.77906 * g / cm3, G4_C, kStateSolid);
0236   }
0237 }
0238
0239 //_______________________________________________________________
0240 void PHG4MicromegasDetector::construct_micromegas(G4LogicalVolume* logicWorld)
0241 {
0242   // start seting up volumes
0243   // Micromegas detector radius
0244   /* it corresponds to the radial position of the innermost surface of a Micromegas module , as measured in CAD model (THREE PANELS UPDATE 1-18-21) */
0245   static constexpr double inner_radius = 83.197 * cm;
0246
0247   /*
0248    * this is the radius at the center of a detector.
0249    * it is updated when constructing the first tile, assuming that all tiles are identical
0250    */
0251   double radius_phi_mean = 0;
0252   double radius_z_mean = 0;
0253
0254   // load survey data
0255   PHG4MicromegasSurvey micromegas_survey;
0256   const bool apply_survey = m_Params->get_int_param("apply_survey");
0257   std::cout << "PHG4MicromegasDetector::construct_micromegas - apply_survey: " << apply_survey << std::endl;
0258
0259   // create detector
0260   // loop over tiles
0261   for (size_t tileid = 0; tileid < m_tiles.size(); ++tileid)
0262   {
0263     // get relevant tile
0264     const auto& tile = m_tiles[tileid];
0265
0266     // create phi tile master volume
0267     auto *const tile_logic_phi = construct_micromegas_tile(tileid, MicromegasDefs::SegmentationType::SEGMENTATION_PHI);
0268     const double tile_thickness_phi = static_cast<G4Box*>(tile_logic_phi->GetSolid())->GetXHalfLength() * 2;
0269
0270     {
0271       // place tile in master volume
0272       const double centerZ = tile.m_centerZ * cm;
0273       const double centerPhi = tile.m_centerPhi;
0274
0275       G4RotationMatrix rotation;
0276       rotation.rotateZ(centerPhi * radian);
0277
0278       // calculate radius at tile center
0279       double radius_phi = inner_radius + tile_thickness_phi / 2;
0280       const G4ThreeVector center(
0281           radius_phi * std::cos(centerPhi),
0282           radius_phi * std::sin(centerPhi),
0283           centerZ);
0284
0285       G4Transform3D transform(rotation, center);
0286
0287       if (apply_survey)
0288       {
0289         // get transformation from survey
0290         const auto survey_transform = micromegas_survey.get_transformation(m_FirstLayer, tileid);
0291         transform = survey_transform * transform;
0292
0293         // adjust radius at tile center to account for survey
0294         const auto translation = transform.getTranslation();
0295         radius_phi = get_r(translation.x(), translation.y());
0296       }
0297
0298       // update mean radius
0299       radius_phi_mean += radius_phi;
0300
0301       const auto tilename = GetName() + "_tile_" + std::to_string(tileid) + "_phi";
0302       new G4PVPlacement(transform, tile_logic_phi, tilename + "_phys", logicWorld, false, 0, OverlapCheck());
0303     }
0304
0305     // create z tile master volume
0306     auto *const tile_logic_z = construct_micromegas_tile(tileid, MicromegasDefs::SegmentationType::SEGMENTATION_Z);
0307     const double tile_thickness_z = static_cast<G4Box*>(tile_logic_z->GetSolid())->GetXHalfLength() * 2;
0308
0309     {
0310       // place tile in master volume
0311       const double centerZ = tile.m_centerZ * cm;
0312       const double centerPhi = tile.m_centerPhi;
0313
0314       G4RotationMatrix rotation;
0315       rotation.rotateZ(centerPhi * radian);
0316
0317       double radius_z = inner_radius + tile_thickness_phi + tile_thickness_z / 2;
0318       const G4ThreeVector center(
0319           radius_z * std::cos(centerPhi),
0320           radius_z * std::sin(centerPhi),
0321           centerZ);
0322
0323       G4Transform3D transform(rotation, center);
0324
0325       if (apply_survey)
0326       {
0327         // get transformation from survey
0328         const auto survey_transform = micromegas_survey.get_transformation(m_FirstLayer + 1, tileid);
0329         transform = survey_transform * transform;
0330
0331         // adjust radius at tile center to account for survey
0332         const auto translation = transform.getTranslation();
0333         radius_z = get_r(translation.x(), translation.y());
0334       }
0335
0336       radius_z_mean += radius_z;
0337
0338       const auto tilename = GetName() + "_tile_" + std::to_string(tileid) + "_z";
0339       new G4PVPlacement(transform, tile_logic_z, tilename + "_phys", logicWorld, false, 0, OverlapCheck());
0340     }
0341   }
0342
0343   // adjust active volume radius to account for world placement
0344   radius_phi_mean /= m_tiles.size();
0345   m_layer_radius.at(m_FirstLayer) += radius_phi_mean / cm;
0346
0347   radius_z_mean /= m_tiles.size();
0348   m_layer_radius.at(m_FirstLayer + 1) += radius_z_mean / cm;
0349
0350   if (Verbosity())
0351   {
0352     std::cout << "PHG4MicromegasDetector::ConstructMe - first layer: " << m_FirstLayer << std::endl;
0353     for (const auto& pair : m_activeVolumes)
0354     {
0355       std::cout << "PHG4MicromegasDetector::ConstructMe - layer: " << pair.second << " volume: " << pair.first->GetName() << std::endl;
0356     }
0357   }
0358
0359   return;
0360 }
0361
0362 //_______________________________________________________________
0363 G4LogicalVolume* PHG4MicromegasDetector::construct_micromegas_tile(int tileid, MicromegasDefs::SegmentationType segmentation_type)
0364 {
0365   // components enumeration
0366   /*
0367   this describes all the detector onion layers for a single side
0368   note that the detector is two sided
0369   */
0370   enum class Component
0371   {
0372     PCB,
0373     CuStrips,
0374     KaptonStrips,
0375     ResistiveStrips,
0376     Gas1,
0377     Mesh,
0378     Gas2,
0379     DriftCuElectrode,
0380     DriftKapton,
0381     DriftCarbon,
0382     FeeSupport
0383   };
0384
0385   // layer definition
0386   struct LayerStruct
0387   {
0388     // constructor
0389     LayerStruct(float thickness, G4Material* material, const G4Colour& color, double dy, double dz, double y_offset, double z_offset)
0390       : m_thickness(thickness)
0391       , m_material(material)
0392       , m_color(color)
0393       , m_dy(dy)
0394       , m_dz(dz)
0395       , m_y_offset(y_offset)
0396       , m_z_offset(z_offset)
0397     {
0398     }
0399 // NOLINTBEGIN(misc-non-private-member-variables-in-classes)
0400     // thickness
0401     float m_thickness = 0;
0402
0403     // material
0404     G4Material* m_material = nullptr;
0405
0406     // color
0407     G4Colour m_color = 0;
0408
0409     // dimension along y
0410     double m_dy = 0;
0411
0412     // dimension along z
0413     double m_dz = 0;
0414
0415     // center offset along y
0416     double m_y_offset = 0;
0417
0418     // center offset along z
0419     double m_z_offset = 0;
0420 // NOLINTEND(misc-non-private-member-variables-in-classes)
0421   };
0422
0423   // define all layers
0424   const std::map<Component, LayerStruct> layer_map =
0425       {
0426           /* adjusted PCB thickness so that total thickness up to Gas2 is 1.6mm, consistently with CAD model */
0427           // { Component::PCB, LayerStruct(1.384*mm, GetDetectorMaterial("mmg_FR4"), G4Colour::Green(), 316*mm, 542*mm, 0, 0 )},
0428           {Component::PCB, LayerStruct(1.0 * mm, GetDetectorMaterial("mmg_FR4"), G4Colour::Green(), 316 * mm, 542 * mm, 0, 0)},
0429           {Component::CuStrips, LayerStruct(12. * micrometer, GetDetectorMaterial("mmg_Strips"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0430           {Component::KaptonStrips, LayerStruct(50. * micrometer, GetDetectorMaterial("mmg_Kapton"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0431           {Component::ResistiveStrips, LayerStruct(20. * micrometer, GetDetectorMaterial("mmg_ResistPaste"), G4Colour::Black(), 256 * mm, 512 * mm, -15 * mm, 0)},
0432           {Component::Gas1, LayerStruct(120. * micrometer, GetDetectorMaterial("mmg_Gas"), G4Colour::Grey(), 256 * mm, 512 * mm, -15 * mm, 0)},
0433           /* 0.8 correction factor to thickness is to account for the mesh denstity@18/45 */
0434           {Component::Mesh, LayerStruct(18. * 0.8 * micrometer, GetDetectorMaterial("mmg_Mesh"), G4Colour::White(), 256 * mm, 512 * mm, -15 * mm, 0)},
0435           {Component::Gas2, LayerStruct(3. * mm, GetDetectorMaterial("mmg_Gas"), G4Colour::Grey(), 256 * mm, 512 * mm, -15 * mm, 0)},
0436           {Component::DriftCuElectrode, LayerStruct(15. * micrometer, GetDetectorMaterial("G4_Cu"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0437           {Component::DriftKapton, LayerStruct(50. * micrometer, GetDetectorMaterial("mmg_Kapton"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0438           {Component::DriftCarbon, LayerStruct(1. * mm, GetDetectorMaterial("G4_C"), G4Colour(150 / 255., 75 / 255., 0), 256 * mm, 512 * mm, -15 * mm, 0)},
0439           {Component::FeeSupport, LayerStruct(2.38 * mm, GetDetectorMaterial("G4_Al"), G4Colour::Grey(), 182 * mm, 542 * mm, -67 * mm, 0)}};
0440
0441   // setup layers in the correct order, going outwards from beam axis
0442   using LayerDefinition = std::tuple<Component, std::string>;
0443   const std::vector<LayerDefinition> layer_stack_phi =
0444       {
0445           // inner side
0446           std::make_tuple(Component::DriftCarbon, "DriftCarbon_inner"),
0447           std::make_tuple(Component::DriftKapton, "DriftKapton_inner"),
0448           std::make_tuple(Component::DriftCuElectrode, "DriftCuElectrode_inner"),
0449           std::make_tuple(Component::Gas2, "Gas2_inner"),
0450           std::make_tuple(Component::Mesh, "Mesh_inner"),
0451           std::make_tuple(Component::Gas1, "Gas1_inner"),
0452           std::make_tuple(Component::ResistiveStrips, "ResistiveStrips_inner"),
0453           std::make_tuple(Component::KaptonStrips, "KaptonStrips_inner"),
0454           std::make_tuple(Component::CuStrips, "CuStrips_inner"),
0455           std::make_tuple(Component::PCB, "PCB_inner")};
0456
0457   // layer stack for z view, same as phi view, but mirrored
0458   const std::vector<LayerDefinition> layer_stack_z =
0459       {
0460           // outer side (= inner side, mirrored)
0461           std::make_tuple(Component::PCB, "PCB_outer"),
0462           std::make_tuple(Component::CuStrips, "CuStrips_outer"),
0463           std::make_tuple(Component::KaptonStrips, "KaptonStrips_outer"),
0464           std::make_tuple(Component::ResistiveStrips, "ResistiveStrips_outer"),
0465           std::make_tuple(Component::Gas1, "Gas1_outer"),
0466           std::make_tuple(Component::Mesh, "Mesh_outer"),
0467           std::make_tuple(Component::Gas2, "Gas2_outer"),
0468           std::make_tuple(Component::DriftCuElectrode, "DriftCuElectrode_outer"),
0469           std::make_tuple(Component::DriftKapton, "DriftKapton_outer"),
0470           std::make_tuple(Component::DriftCarbon, "DriftCarbon_outer"),
0471
0472           // FEE support plate
0473           std::make_tuple(Component::FeeSupport, "FEE_support")};
0474
0475   const bool is_z = (segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_Z);
0476
0477   // get the right layer stack
0478   const auto& layer_stack = is_z ? layer_stack_z : layer_stack_phi;
0479
0480   // for z view, create two FEE boards up front, to get their total thickness and add to master volume
0481   std::array<G4LogicalVolume*, 2> fee_board_logic =
0482       {
0483           is_z ? construct_fee_board(0) : nullptr,
0484           is_z ? construct_fee_board(1) : nullptr};
0485
0486   const double fee_thickness = is_z ? static_cast<G4Box*>(fee_board_logic[0]->GetSolid())->GetXHalfLength() * 2 : 0;
0487
0488   // calculate total tile thickness
0489   double tile_thickness = std::accumulate(
0490       layer_stack.begin(), layer_stack.end(), 0.,
0491       [&layer_map](double value, const LayerDefinition& layer)
0492       { return value + layer_map.at(std::get<0>(layer)).m_thickness; });
0493
0494   // for z tile, adds fee thickness
0495   if (is_z)
0496   {
0497     tile_thickness += fee_thickness;
0498   }
0499
0500   // tile dimensions match that of the PCB layer
0501   const double tile_dy = layer_map.at(Component::PCB).m_dy;
0502   const double tile_dz = layer_map.at(Component::PCB).m_dz;
0503
0504   // get world material to define parent volume
0505   auto *rc = recoConsts::instance();
0506   auto *world_material = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0507
0508   // define tile name
0509   const auto tilename = GetName() + "_tile_" + std::to_string(tileid) + (is_z ? "_z" : "_phi");
0510
0511   auto *tile_solid = new G4Box(tilename + "_solid", tile_thickness / 2, tile_dy / 2, tile_dz / 2);
0512   auto *tile_logic = new G4LogicalVolume(tile_solid, world_material, "invisible_" + tilename + "_logic");
0513   GetDisplayAction()->AddVolume(tile_logic, G4Colour::Grey());
0514
0515   /* we loop over registered layers and create volumes for each as daughter of the tile volume */
0516   auto current_radius_local = -tile_thickness / 2;
0517   for (const auto& [type, name] : layer_stack)
0518   {
0519     // layer name
0520     /*
0521      * for the Gas2 layers, which are the active components, we use a different volume name,
0522      * that match the old geometry implementation. This maximizes compatibility with previous versions
0523      */
0524     const G4String cname = (type == Component::Gas2) ? "micromegas_measurement_" + name : G4String(GetName()) + "_" + name;
0525
0526     // get thickness, material and name
0527     const auto& component(layer_map.at(type));
0528     const auto& thickness = component.m_thickness;
0529     const auto& material = component.m_material;
0530     const auto& color = component.m_color;
0531
0532     const auto& dy = component.m_dy;
0533     const auto& dz = component.m_dz;
0534
0535     const auto& y_offset = component.m_y_offset;
0536     const auto& z_offset = component.m_z_offset;
0537
0538     auto *component_solid = new G4Box(cname + "_solid", thickness / 2, dy / 2, dz / 2);
0539     auto *component_logic = new G4LogicalVolume(component_solid, material, cname + "_logic");
0540     GetDisplayAction()->AddVolume(component_logic, color);
0541
0542     const G4ThreeVector center((current_radius_local + thickness / 2), y_offset, z_offset);
0543     auto *component_phys = new G4PVPlacement(nullptr, center, component_logic, cname + "_phys", tile_logic, false, 0, OverlapCheck());
0544
0545     if (type == Component::Gas2)
0546     {
0547       // store active volume
0548       // define layer from name
0549       const int layer_index = is_z ? m_FirstLayer + 1 : m_FirstLayer;
0550       m_activeVolumes.insert(std::make_pair(component_phys, layer_index));
0551       m_tiles_map.insert(std::make_pair(component_phys, tileid));
0552
0553       // store radius associated to this layer
0554       m_layer_radius.insert(std::make_pair(layer_index, (current_radius_local + thickness / 2) / cm));
0555       m_layer_thickness.insert(std::make_pair(layer_index, thickness / cm));
0556     }
0557     else
0558     {
0559       m_passiveVolumes.insert(component_phys);
0560     }
0561
0562     // update radius
0563     current_radius_local += thickness;
0564   }
0565
0566   if (is_z)
0567   {
0568     // add FEE boards
0569     /* offsets measured from CAD drawings */
0570     static constexpr double fee_y_offset = (316. / 2 - 44.7 - 141.5 / 2) * mm;
0571     static constexpr double fee_x_offset = (542. / 2 - 113.1 - 140. / 2) * mm;
0572     new G4PVPlacement(nullptr, {current_radius_local + fee_thickness / 2, -fee_y_offset, -fee_x_offset}, fee_board_logic[0], GetName() + "_fee_0_phys", tile_logic, false, 0, OverlapCheck());
0573     new G4PVPlacement(nullptr, {current_radius_local + fee_thickness / 2, -fee_y_offset, fee_x_offset}, fee_board_logic[1], GetName() + "_fee_1_phys", tile_logic, false, 0, OverlapCheck());
0574   }
0575
0576   // return master logical volume
0577   return tile_logic;
0578 }
0579
0580 //_______________________________________________________________
0581 G4LogicalVolume* PHG4MicromegasDetector::construct_fee_board(int id)
0582 {
0583   // layer definition
0584   struct LayerStruct
0585   {
0586     // constructor
0587     LayerStruct(const std::string& name, float thickness, G4Material* material, const G4Colour& color)
0588       : m_name(name)
0589       , m_thickness(thickness)
0590       , m_material(material)
0591       , m_color(color)
0592     {
0593     }
0594
0595 // NOLINTBEGIN(misc-non-private-member-variables-in-classes)
0596     // name
0597     std::string m_name;
0598
0599     // thickness
0600     float m_thickness = 0;
0601
0602     // material
0603     G4Material* m_material = nullptr;
0604
0605     // color
0606     G4Colour m_color = 0;
0607 // NOLINTEND(misc-non-private-member-variables-in-classes)
0608   };
0609
0610   static constexpr double inch_to_cm = 2.54;
0611
0612   /*
0613    * FEE board consists of FR4 PCB, a coper layer, and an aluminium layer, for cooling.
0614    * FR4 and Cu layer thickness taken from TPC description (PHG4TpcEndCapSubsystem::SetDefaultParameters)
0615    * Al layer correspond to cooling plate. Thickness from mechanical drawings
0616    */
0617   const std::vector<LayerStruct> layer_stack = {
0618       LayerStruct("fee_pcb", 0.07 * inch_to_cm * cm, GetDetectorMaterial("mmg_FR4"), G4Colour::Green()),
0619       LayerStruct("fee_cu", 35e-4 * 10 * 0.8 * cm, GetDetectorMaterial("G4_Cu"), G4Colour::Brown()),
0620       LayerStruct("fee_al", 0.25 * inch_to_cm * cm, GetDetectorMaterial("G4_Al"), G4Colour::Grey())};
0621
0622   // calculate total tile thickness
0623   const double fee_thickness = std::accumulate(
0624       layer_stack.begin(), layer_stack.end(), 0.,
0625       [](double value, const LayerStruct& layer)
0626       { return value + layer.m_thickness; });
0627
0628   // fee dimensions match CAD drawings
0629   const double fee_dy = 141.51 * mm;
0630   const double fee_dz = 140 * mm;
0631
0632   // get world material to define parent volume
0633   auto *rc = recoConsts::instance();
0634   auto *world_material = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0635
0636   // define tile name
0637   const auto feename = GetName() + "_fee_" + std::to_string(id);
0638
0639   auto *fee_solid = new G4Box(feename + "_solid", fee_thickness / 2, fee_dy / 2, fee_dz / 2);
0640   auto *fee_logic = new G4LogicalVolume(fee_solid, world_material, "invisible_" + feename + "_logic");
0641   GetDisplayAction()->AddVolume(fee_logic, G4Colour::Grey());
0642
0643   /* we loop over registered layers and create volumes for each as daughter of the fee volume */
0644   auto current_radius_local = -fee_thickness / 2;
0645   for (const auto& layer : layer_stack)
0646   {
0647     // layer name
0648     /*
0649      * for the Gas2 layers, which are the active components, we use a different volume name,
0650      * that match the old geometry implementation. This maximizes compatibility with previous versions
0651      */
0652     const G4String cname = G4String(GetName()) + "_" + layer.m_name;
0653
0654     // get thickness, material and name
0655     const auto& thickness = layer.m_thickness;
0656     const auto& material = layer.m_material;
0657     const auto& color = layer.m_color;
0658
0659     auto *component_solid = new G4Box(cname + "_solid", thickness / 2, fee_dy / 2, fee_dz / 2);
0660     auto *component_logic = new G4LogicalVolume(component_solid, material, cname + "_logic");
0661     GetDisplayAction()->AddVolume(component_logic, color);
0662
0663     const G4ThreeVector center((current_radius_local + thickness / 2), 0, 0);
0664     auto *component_phys = new G4PVPlacement(nullptr, center, component_logic, cname + "_phys", fee_logic, false, 0, OverlapCheck());
0665
0666     // store as passive
0667     m_passiveVolumes.insert(component_phys);
0668
0669     // update radius
0670     current_radius_local += thickness;
0671   }
0672
0673   // return master logical volume
0674   return fee_logic;
0675 }
0676
0677 //_______________________________________________________________
0678 void PHG4MicromegasDetector::add_geometry_node()
0679 {
0680   // do nothing if detector is inactive
0681   if (!m_Params->get_int_param("active"))
0682   {
0683     return;
0684   }
0685
0686   // find or create geometry node
0687   const auto geonode_name = std::string("CYLINDERGEOM_") + m_SuperDetector + "_FULL";
0688   auto *geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode_name);
0689   if (!geonode)
0690   {
0691     geonode = new PHG4CylinderGeomContainer;
0692     PHNodeIterator iter(topNode());
0693     auto *runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0694     PHNodeIterator runiter(runNode);
0695     PHCompositeNode *geomNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", "RECO_TRACKING_GEOMETRY"));
0696     if(!geomNode)
0697       {
0698         geomNode = new PHCompositeNode("RECO_TRACKING_GEOMETRY");
0699         runNode->addNode(geomNode);
0700       }
0701     auto *newNode = new PHIODataNode<PHObject>(geonode, geonode_name, "PHObject");
0702     geomNode->addNode(newNode);
0703   }
0704
0705   // cylinder maximal length
0706   static constexpr double length = 220;
0707
0708   // add cylinder objects
0709   /* one cylinder is added per physical volume. The dimention correspond to the drift volume */
0710   bool is_first = true;
0711   for (const auto& [layer_index, radius] : m_layer_radius)
0712   {
0713     // create cylinder and match geometry
0714     /* note: cylinder segmentation type and pitch is set in PHG4MicromegasHitReco */
0715     auto *cylinder = new CylinderGeomMicromegas(layer_index);
0716     cylinder->set_radius(radius);
0717     cylinder->set_thickness(m_layer_thickness.at(layer_index));
0718     cylinder->set_zmin(-length / 2);
0719     cylinder->set_zmax(length / 2);
0720
0721     // tiles
0722     cylinder->set_tiles(m_tiles);
0723
0724     /*
0725      * asign segmentation type and pitch
0726      * assume first layer in phi, other(s) are z
0727      */
0728     cylinder->set_segmentation_type(is_first ? MicromegasDefs::SegmentationType::SEGMENTATION_PHI : MicromegasDefs::SegmentationType::SEGMENTATION_Z);
0729
0730     /*
0731      * assign drift direction
0732      * assume first layer is outward, with readout plane at the top, and others are inward, with readout plane at the bottom
0733      * this is used to properly implement transverse diffusion in ::process_event
0734      */
0735     cylinder->set_drift_direction(is_first ? MicromegasDefs::DriftDirection::OUTWARD : MicromegasDefs::DriftDirection::INWARD);
0736
0737     // pitch
0738     // first layer (phi) has 1mm pitch. second layer (z) has 2mm pitch
0739     cylinder->set_pitch(is_first ? 0.1 : 0.2);
0740
0741     // if( Verbosity() )
0742     {
0743       cylinder->identify(std::cout);
0744     }
0745
0746     geonode->AddLayerGeom(layer_index, cylinder);
0747
0748     is_first = false;
0749   }
0750 }