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   inline 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.find(volume) != m_activeVolumes.end())
0080     {
0081       return 1;
0082     }
0083   }
0084   if (m_SupportActiveFlag)
0085   {
0086     if (m_passiveVolumes.find(volume) != m_passiveVolumes.end())
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     const auto 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     const auto 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
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   };
0421
0422   // define all layers
0423   const std::map<Component, LayerStruct> layer_map =
0424       {
0425           /* adjusted PCB thickness so that total thickness up to Gas2 is 1.6mm, consistently with CAD model */
0426           // { Component::PCB, LayerStruct(1.384*mm, GetDetectorMaterial("mmg_FR4"), G4Colour::Green(), 316*mm, 542*mm, 0, 0 )},
0427           {Component::PCB, LayerStruct(1.0 * mm, GetDetectorMaterial("mmg_FR4"), G4Colour::Green(), 316 * mm, 542 * mm, 0, 0)},
0428           {Component::CuStrips, LayerStruct(12. * micrometer, GetDetectorMaterial("mmg_Strips"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0429           {Component::KaptonStrips, LayerStruct(50. * micrometer, GetDetectorMaterial("mmg_Kapton"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0430           {Component::ResistiveStrips, LayerStruct(20. * micrometer, GetDetectorMaterial("mmg_ResistPaste"), G4Colour::Black(), 256 * mm, 512 * mm, -15 * mm, 0)},
0431           {Component::Gas1, LayerStruct(120. * micrometer, GetDetectorMaterial("mmg_Gas"), G4Colour::Grey(), 256 * mm, 512 * mm, -15 * mm, 0)},
0432           /* 0.8 correction factor to thickness is to account for the mesh denstity@18/45 */
0433           {Component::Mesh, LayerStruct(18. * 0.8 * micrometer, GetDetectorMaterial("mmg_Mesh"), G4Colour::White(), 256 * mm, 512 * mm, -15 * mm, 0)},
0434           {Component::Gas2, LayerStruct(3. * mm, GetDetectorMaterial("mmg_Gas"), G4Colour::Grey(), 256 * mm, 512 * mm, -15 * mm, 0)},
0435           {Component::DriftCuElectrode, LayerStruct(15. * micrometer, GetDetectorMaterial("G4_Cu"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0436           {Component::DriftKapton, LayerStruct(50. * micrometer, GetDetectorMaterial("mmg_Kapton"), G4Colour::Brown(), 256 * mm, 512 * mm, -15 * mm, 0)},
0437           {Component::DriftCarbon, LayerStruct(1. * mm, GetDetectorMaterial("G4_C"), G4Colour(150 / 255., 75 / 255., 0), 256 * mm, 512 * mm, -15 * mm, 0)},
0438           {Component::FeeSupport, LayerStruct(2.38 * mm, GetDetectorMaterial("G4_Al"), G4Colour::Grey(), 182 * mm, 542 * mm, -67 * mm, 0)}};
0439
0440   // setup layers in the correct order, going outwards from beam axis
0441   using LayerDefinition = std::tuple<Component, std::string>;
0442   const std::vector<LayerDefinition> layer_stack_phi =
0443       {
0444           // inner side
0445           std::make_tuple(Component::DriftCarbon, "DriftCarbon_inner"),
0446           std::make_tuple(Component::DriftKapton, "DriftKapton_inner"),
0447           std::make_tuple(Component::DriftCuElectrode, "DriftCuElectrode_inner"),
0448           std::make_tuple(Component::Gas2, "Gas2_inner"),
0449           std::make_tuple(Component::Mesh, "Mesh_inner"),
0450           std::make_tuple(Component::Gas1, "Gas1_inner"),
0451           std::make_tuple(Component::ResistiveStrips, "ResistiveStrips_inner"),
0452           std::make_tuple(Component::KaptonStrips, "KaptonStrips_inner"),
0453           std::make_tuple(Component::CuStrips, "CuStrips_inner"),
0454           std::make_tuple(Component::PCB, "PCB_inner")};
0455
0456   // layer stack for z view, same as phi view, but mirrored
0457   const std::vector<LayerDefinition> layer_stack_z =
0458       {
0459           // outer side (= inner side, mirrored)
0460           std::make_tuple(Component::PCB, "PCB_outer"),
0461           std::make_tuple(Component::CuStrips, "CuStrips_outer"),
0462           std::make_tuple(Component::KaptonStrips, "KaptonStrips_outer"),
0463           std::make_tuple(Component::ResistiveStrips, "ResistiveStrips_outer"),
0464           std::make_tuple(Component::Gas1, "Gas1_outer"),
0465           std::make_tuple(Component::Mesh, "Mesh_outer"),
0466           std::make_tuple(Component::Gas2, "Gas2_outer"),
0467           std::make_tuple(Component::DriftCuElectrode, "DriftCuElectrode_outer"),
0468           std::make_tuple(Component::DriftKapton, "DriftKapton_outer"),
0469           std::make_tuple(Component::DriftCarbon, "DriftCarbon_outer"),
0470
0471           // FEE support plate
0472           std::make_tuple(Component::FeeSupport, "FEE_support")};
0473
0474   const bool is_z = (segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_Z);
0475
0476   // get the right layer stack
0477   const auto& layer_stack = is_z ? layer_stack_z : layer_stack_phi;
0478
0479   // for z view, create two FEE boards up front, to get their total thickness and add to master volume
0480   std::array<G4LogicalVolume*, 2> fee_board_logic =
0481       {
0482           is_z ? construct_fee_board(0) : nullptr,
0483           is_z ? construct_fee_board(1) : nullptr};
0484
0485   const double fee_thickness = is_z ? static_cast<G4Box*>(fee_board_logic[0]->GetSolid())->GetXHalfLength() * 2 : 0;
0486
0487   // calculate total tile thickness
0488   double tile_thickness = std::accumulate(
0489       layer_stack.begin(), layer_stack.end(), 0.,
0490       [&layer_map](double value, const LayerDefinition& layer)
0491       { return value + layer_map.at(std::get<0>(layer)).m_thickness; });
0492
0493   // for z tile, adds fee thickness
0494   if (is_z)
0495   {
0496     tile_thickness += fee_thickness;
0497   }
0498
0499   // tile dimensions match that of the PCB layer
0500   const double tile_dy = layer_map.at(Component::PCB).m_dy;
0501   const double tile_dz = layer_map.at(Component::PCB).m_dz;
0502
0503   // get world material to define parent volume
0504   auto rc = recoConsts::instance();
0505   auto world_material = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0506
0507   // define tile name
0508   const auto tilename = GetName() + "_tile_" + std::to_string(tileid) + (is_z ? "_z" : "_phi");
0509
0510   auto tile_solid = new G4Box(tilename + "_solid", tile_thickness / 2, tile_dy / 2, tile_dz / 2);
0511   auto tile_logic = new G4LogicalVolume(tile_solid, world_material, "invisible_" + tilename + "_logic");
0512   GetDisplayAction()->AddVolume(tile_logic, G4Colour::Grey());
0513
0514   /* we loop over registered layers and create volumes for each as daughter of the tile volume */
0515   auto current_radius_local = -tile_thickness / 2;
0516   for (const auto& [type, name] : layer_stack)
0517   {
0518     // layer name
0519     /*
0520      * for the Gas2 layers, which are the active components, we use a different volume name,
0521      * that match the old geometry implementation. This maximizes compatibility with previous versions
0522      */
0523     const G4String cname = (type == Component::Gas2) ? "micromegas_measurement_" + name : G4String(GetName()) + "_" + name;
0524
0525     // get thickness, material and name
0526     const auto& component(layer_map.at(type));
0527     const auto& thickness = component.m_thickness;
0528     const auto& material = component.m_material;
0529     const auto& color = component.m_color;
0530
0531     const auto& dy = component.m_dy;
0532     const auto& dz = component.m_dz;
0533
0534     const auto& y_offset = component.m_y_offset;
0535     const auto& z_offset = component.m_z_offset;
0536
0537     auto component_solid = new G4Box(cname + "_solid", thickness / 2, dy / 2, dz / 2);
0538     auto component_logic = new G4LogicalVolume(component_solid, material, cname + "_logic");
0539     GetDisplayAction()->AddVolume(component_logic, color);
0540
0541     const G4ThreeVector center((current_radius_local + thickness / 2), y_offset, z_offset);
0542     auto component_phys = new G4PVPlacement(nullptr, center, component_logic, cname + "_phys", tile_logic, false, 0, OverlapCheck());
0543
0544     if (type == Component::Gas2)
0545     {
0546       // store active volume
0547       // define layer from name
0548       const int layer_index = is_z ? m_FirstLayer + 1 : m_FirstLayer;
0549       m_activeVolumes.insert(std::make_pair(component_phys, layer_index));
0550       m_tiles_map.insert(std::make_pair(component_phys, tileid));
0551
0552       // store radius associated to this layer
0553       m_layer_radius.insert(std::make_pair(layer_index, (current_radius_local + thickness / 2) / cm));
0554       m_layer_thickness.insert(std::make_pair(layer_index, thickness / cm));
0555     }
0556     else
0557     {
0558       m_passiveVolumes.insert(component_phys);
0559     }
0560
0561     // update radius
0562     current_radius_local += thickness;
0563   }
0564
0565   if (is_z)
0566   {
0567     // add FEE boards
0568     /* offsets measured from CAD drawings */
0569     static constexpr double fee_y_offset = (316. / 2 - 44.7 - 141.5 / 2) * mm;
0570     static constexpr double fee_x_offset = (542. / 2 - 113.1 - 140. / 2) * mm;
0571     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());
0572     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());
0573   }
0574
0575   // return master logical volume
0576   return tile_logic;
0577 }
0578
0579 //_______________________________________________________________
0580 G4LogicalVolume* PHG4MicromegasDetector::construct_fee_board(int id)
0581 {
0582   // layer definition
0583   struct LayerStruct
0584   {
0585     // constructor
0586     LayerStruct(const std::string& name, float thickness, G4Material* material, const G4Colour& color)
0587       : m_name(name)
0588       , m_thickness(thickness)
0589       , m_material(material)
0590       , m_color(color)
0591     {
0592     }
0593
0594     // name
0595     std::string m_name;
0596
0597     // thickness
0598     float m_thickness = 0;
0599
0600     // material
0601     G4Material* m_material = nullptr;
0602
0603     // color
0604     G4Colour m_color = 0;
0605   };
0606
0607   static constexpr double inch_to_cm = 2.54;
0608
0609   /*
0610    * FEE board consists of FR4 PCB, a coper layer, and an aluminium layer, for cooling.
0611    * FR4 and Cu layer thickness taken from TPC description (PHG4TpcEndCapSubsystem::SetDefaultParameters)
0612    * Al layer correspond to cooling plate. Thickness from mechanical drawings
0613    */
0614   const std::vector<LayerStruct> layer_stack = {
0615       LayerStruct("fee_pcb", 0.07 * inch_to_cm * cm, GetDetectorMaterial("mmg_FR4"), G4Colour::Green()),
0616       LayerStruct("fee_cu", 35e-4 * 10 * 0.8 * cm, GetDetectorMaterial("G4_Cu"), G4Colour::Brown()),
0617       LayerStruct("fee_al", 0.25 * inch_to_cm * cm, GetDetectorMaterial("G4_Al"), G4Colour::Grey())};
0618
0619   // calculate total tile thickness
0620   const double fee_thickness = std::accumulate(
0621       layer_stack.begin(), layer_stack.end(), 0.,
0622       [](double value, const LayerStruct& layer)
0623       { return value + layer.m_thickness; });
0624
0625   // fee dimensions match CAD drawings
0626   const double fee_dy = 141.51 * mm;
0627   const double fee_dz = 140 * mm;
0628
0629   // get world material to define parent volume
0630   auto rc = recoConsts::instance();
0631   auto world_material = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0632
0633   // define tile name
0634   const auto feename = GetName() + "_fee_" + std::to_string(id);
0635
0636   auto fee_solid = new G4Box(feename + "_solid", fee_thickness / 2, fee_dy / 2, fee_dz / 2);
0637   auto fee_logic = new G4LogicalVolume(fee_solid, world_material, "invisible_" + feename + "_logic");
0638   GetDisplayAction()->AddVolume(fee_logic, G4Colour::Grey());
0639
0640   /* we loop over registered layers and create volumes for each as daughter of the fee volume */
0641   auto current_radius_local = -fee_thickness / 2;
0642   for (const auto& layer : layer_stack)
0643   {
0644     // layer name
0645     /*
0646      * for the Gas2 layers, which are the active components, we use a different volume name,
0647      * that match the old geometry implementation. This maximizes compatibility with previous versions
0648      */
0649     const G4String cname = G4String(GetName()) + "_" + layer.m_name;
0650
0651     // get thickness, material and name
0652     const auto& thickness = layer.m_thickness;
0653     const auto& material = layer.m_material;
0654     const auto& color = layer.m_color;
0655
0656     auto component_solid = new G4Box(cname + "_solid", thickness / 2, fee_dy / 2, fee_dz / 2);
0657     auto component_logic = new G4LogicalVolume(component_solid, material, cname + "_logic");
0658     GetDisplayAction()->AddVolume(component_logic, color);
0659
0660     const G4ThreeVector center((current_radius_local + thickness / 2), 0, 0);
0661     auto component_phys = new G4PVPlacement(nullptr, center, component_logic, cname + "_phys", fee_logic, false, 0, OverlapCheck());
0662
0663     // store as passive
0664     m_passiveVolumes.insert(component_phys);
0665
0666     // update radius
0667     current_radius_local += thickness;
0668   }
0669
0670   // return master logical volume
0671   return fee_logic;
0672 }
0673
0674 //_______________________________________________________________
0675 void PHG4MicromegasDetector::add_geometry_node()
0676 {
0677   // do nothing if detector is inactive
0678   if (!m_Params->get_int_param("active"))
0679   {
0680     return;
0681   }
0682
0683   // find or create geometry node
0684   const auto geonode_name = std::string("CYLINDERGEOM_") + m_SuperDetector + "_FULL";
0685   auto geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode_name);
0686   if (!geonode)
0687   {
0688     geonode = new PHG4CylinderGeomContainer;
0689     PHNodeIterator iter(topNode());
0690     auto runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0691     auto newNode = new PHIODataNode<PHObject>(geonode, geonode_name, "PHObject");
0692     runNode->addNode(newNode);
0693   }
0694
0695   // cylinder maximal length
0696   static constexpr double length = 220;
0697
0698   // add cylinder objects
0699   /* one cylinder is added per physical volume. The dimention correspond to the drift volume */
0700   bool is_first = true;
0701   for (const auto& [layer_index, radius] : m_layer_radius)
0702   {
0703     // create cylinder and match geometry
0704     /* note: cylinder segmentation type and pitch is set in PHG4MicromegasHitReco */
0705     auto cylinder = new CylinderGeomMicromegas(layer_index);
0706     cylinder->set_radius(radius);
0707     cylinder->set_thickness(m_layer_thickness.at(layer_index));
0708     cylinder->set_zmin(-length / 2);
0709     cylinder->set_zmax(length / 2);
0710
0711     // tiles
0712     cylinder->set_tiles(m_tiles);
0713
0714     /*
0715      * asign segmentation type and pitch
0716      * assume first layer in phi, other(s) are z
0717      */
0718     cylinder->set_segmentation_type(is_first ? MicromegasDefs::SegmentationType::SEGMENTATION_PHI : MicromegasDefs::SegmentationType::SEGMENTATION_Z);
0719
0720     /*
0721      * assign drift direction
0722      * assume first layer is outward, with readout plane at the top, and others are inward, with readout plane at the bottom
0723      * this is used to properly implement transverse diffusion in ::process_event
0724      */
0725     cylinder->set_drift_direction(is_first ? MicromegasDefs::DriftDirection::OUTWARD : MicromegasDefs::DriftDirection::INWARD);
0726
0727     // pitch
0728     // first layer (phi) has 1mm pitch. second layer (z) has 2mm pitch
0729     cylinder->set_pitch(is_first ? 0.1 : 0.2);
0730
0731     // if( Verbosity() )
0732     {
0733       cylinder->identify(std::cout);
0734     }
0735
0736     geonode->AddLayerGeom(layer_index, cylinder);
0737
0738     is_first = false;
0739   }
0740 }