g4simulation/g4detectors/PHG4SectorConstructor.cc

0001 /*!
0002  * \file PHG4SectorConstructor.cc
0003  * \brief
0004  * \author Jin Huang <jhuang@bnl.gov>
0005  * \version $Revision: 1.6 $
0006  * \date $Date: 2014/07/31 15:52:37 $
0007  */
0008
0009 #include "PHG4SectorConstructor.h"
0010 #include "PHG4SectorDisplayAction.h"
0011
0012 #include <g4main/PHG4Detector.h>
0013 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0014 #include <g4main/PHG4Subsystem.h>      // for PHG4Subsystem
0015
0016 #include <Geant4/G4Box.hh>
0017 #include <Geant4/G4DisplacedSolid.hh>     // for G4DisplacedSolid
0018 #include <Geant4/G4Exception.hh>          // for G4Exception
0019 #include <Geant4/G4ExceptionSeverity.hh>  // for FatalException, JustWarning
0020 #include <Geant4/G4IntersectionSolid.hh>
0021 #include <Geant4/G4LogicalVolume.hh>
0022 #include <Geant4/G4PVPlacement.hh>
0023 #include <Geant4/G4PhysicalConstants.hh>  // for pi
0024 #include <Geant4/G4Sphere.hh>
0025 #include <Geant4/G4String.hh>
0026 #include <Geant4/G4SystemOfUnits.hh>  // for cm, um, perCent
0027 #include <Geant4/G4ThreeVector.hh>    // for G4ThreeVector
0028 #include <Geant4/G4Transform3D.hh>    // for G4Transform3D, G4RotateX3D
0029 #include <Geant4/G4Tubs.hh>
0030 #include <Geant4/G4Types.hh>  // for G4int
0031
0032 #include <algorithm>  // for max
0033 #include <cassert>
0034 #include <climits>
0035 #include <cmath>
0036 #include <iostream>
0037 #include <sstream>
0038
0039 class G4Material;
0040
0041 PHG4Sector::PHG4SectorConstructor::PHG4SectorConstructor(const std::string &name, PHG4Subsystem *subsys)
0042   : overlapcheck_sector(false)
0043   , name_base(name)
0044   , m_DisplayAction(dynamic_cast<PHG4SectorDisplayAction *>(subsys->GetDisplayAction()))
0045   , m_Verbosity(0)
0046 {
0047 }
0048
0049 void PHG4Sector::PHG4SectorConstructor::Construct_Sectors(G4LogicalVolume *WorldLog)
0050 {
0051   // geometry checks
0052   if (geom.get_total_thickness() == 0)
0053   {
0054     G4Exception(
0055         (std::string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
0056         __FILE__, FatalException,
0057         " detector configured with zero thickness!");
0058   }
0059
0060   if (geom.get_min_polar_angle() == geom.get_max_polar_angle())
0061   {
0062     G4Exception(
0063         (std::string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
0064         __FILE__, FatalException, "min_polar_angle = max_polar_angle!");
0065   }
0066   if (geom.get_min_polar_angle() > geom.get_max_polar_angle())
0067   {
0068     G4Exception(
0069         (std::string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
0070         __FILE__, JustWarning,
0071         "min and max polar angle got reversed. Correcting them.");
0072
0073     const double t = geom.get_max_polar_angle();
0074     geom.set_max_polar_angle(geom.get_min_polar_angle());
0075     geom.set_min_polar_angle(t);
0076   }
0077   if ((geom.get_min_polar_edge() == Sector_Geometry::kFlatEdge or geom.get_max_polar_edge() == Sector_Geometry::kFlatEdge) and geom.get_N_Sector() <= 2)
0078   {
0079     G4Exception(
0080         (std::string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
0081         __FILE__, FatalException,
0082         "can NOT use flat edge for single or double sector detector!");
0083   }
0084
0085   const G4Transform3D transform_Det_to_Hall =
0086       G4RotateX3D(-geom.get_normal_polar_angle()) * G4TranslateZ3D(
0087                                                         geom.get_normal_start() + geom.get_total_thickness() / 2);
0088
0089   const G4Transform3D transform_Hall_to_Det(transform_Det_to_Hall.inverse());
0090
0091   // during GDML export, numerical value may change at the large digit and go beyond the 0 or pi limit.
0092   // therefore recess 0/pi by a small amount to avoid such problem
0093   static const double epsilon = std::numeric_limits<float>::epsilon();
0094   const double sph_min_polar_angle =
0095       (geom.get_min_polar_edge() == Sector_Geometry::kConeEdge) ? geom.get_min_polar_angle() : (0 + epsilon);
0096   const double sph_max_polar_angle =
0097       (geom.get_max_polar_edge() == Sector_Geometry::kConeEdge) ? geom.get_max_polar_angle() : (pi - epsilon);
0098
0099   G4VSolid *SecConeBoundary_Hall = new G4Sphere("SecConeBoundary_Hall",                                           //
0100                                                 geom.get_normal_start(), geom.get_max_R(),                        // G4double pRmin, G4double pRmax,
0101                                                 pi / 2 - pi / geom.get_N_Sector(), 2 * pi / geom.get_N_Sector(),  //  G4double pSPhi, G4double pDPhi,
0102                                                 sph_min_polar_angle, sph_max_polar_angle - sph_min_polar_angle    // G4double pSTheta, G4double pDTheta
0103   );
0104
0105   G4VSolid *SecConeBoundary_Det = new G4DisplacedSolid("SecConeBoundary_Det",
0106                                                        SecConeBoundary_Hall, transform_Hall_to_Det);
0107
0108   G4VSolid *Boundary_Det = SecConeBoundary_Det;
0109
0110   if (geom.get_min_polar_edge() != Sector_Geometry::kConeEdge or geom.get_max_polar_edge() != Sector_Geometry::kConeEdge)
0111   {
0112     // build a flat edge
0113
0114     const double sph_min_polar_size =
0115         (geom.get_min_polar_edge() == Sector_Geometry::kConeEdge) ? geom.get_max_R() : geom.get_normal_start() * tan(geom.get_normal_polar_angle() - geom.get_min_polar_angle());
0116     const double sph_max_polar_size =
0117         (geom.get_max_polar_edge() == Sector_Geometry::kConeEdge) ? geom.get_max_R() : geom.get_normal_start() * tan(geom.get_max_polar_angle() - geom.get_normal_polar_angle());
0118
0119     G4VSolid *BoxBoundary_Det = new G4Box("BoxBoundary_Det",
0120                                           geom.get_max_R(), (sph_min_polar_size + sph_max_polar_size) / 2,
0121                                           geom.get_total_thickness());
0122     G4VSolid *BoxBoundary_Det_Place = new G4DisplacedSolid(
0123         "BoxBoundary_Det_Place", BoxBoundary_Det, nullptr,
0124         G4ThreeVector(0, (sph_max_polar_size - sph_min_polar_size) / 2, 0));
0125
0126     Boundary_Det = new G4IntersectionSolid("Boundary_Det",
0127                                            BoxBoundary_Det_Place, SecConeBoundary_Det);
0128   }
0129
0130   G4VSolid *DetectorBox_Det = Construct_Sectors_Plane("DetectorBox_Det",
0131                                                       -geom.get_total_thickness() / 2, geom.get_total_thickness(),
0132                                                       Boundary_Det);
0133
0134   G4Material *p_mat = PHG4Detector::GetDetectorMaterial(geom.get_material());
0135
0136   G4LogicalVolume *DetectorLog_Det = new G4LogicalVolume(DetectorBox_Det,  //
0137                                                          p_mat, name_base + "_Log");
0138   RegisterLogicalVolume(DetectorLog_Det);
0139
0140   for (G4int sec = 0; sec < geom.get_N_Sector(); sec++)
0141   {
0142     RegisterPhysicalVolume(
0143         new G4PVPlacement(
0144             G4RotateZ3D(2 * pi / geom.get_N_Sector() * sec) * transform_Det_to_Hall, DetectorLog_Det,
0145             name_base + "_Physical", WorldLog, false, sec, overlapcheck_sector));
0146   }
0147
0148   // construct layers
0149   double z_start = -geom.get_total_thickness() / 2;
0150
0151   for (const auto &l : geom.layer_list)
0152   {
0153     if (l.percentage_filled > 100. || l.percentage_filled < 0)
0154     {
0155       std::ostringstream strstr;
0156       strstr << name_base << " have invalid layer ";
0157       strstr << l.name << " with percentage_filled =" << l.percentage_filled;
0158
0159       G4Exception(
0160           (std::string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(), __FILE__, FatalException,
0161           strstr.str().c_str());
0162     }
0163
0164     const std::string layer_name = name_base + "_" + l.name;
0165
0166     G4VSolid *LayerSol_Det = Construct_Sectors_Plane(layer_name + "_Sol",
0167                                                      z_start, l.depth * l.percentage_filled * perCent, Boundary_Det);
0168
0169     G4LogicalVolume *LayerLog_Det = new G4LogicalVolume(LayerSol_Det,  //
0170                                                         PHG4Detector::GetDetectorMaterial(l.material), layer_name + "_Log");
0171     RegisterLogicalVolume(LayerLog_Det);
0172
0173     RegisterPhysicalVolume(
0174         new G4PVPlacement(nullptr, G4ThreeVector(), LayerLog_Det,
0175                           layer_name + "_Physical", DetectorLog_Det, false, 0, overlapcheck_sector),
0176         l.active);
0177
0178     z_start += l.depth;
0179   }
0180
0181   if (std::abs(z_start - geom.get_total_thickness() / 2) > 1 * um)
0182   {
0183     std::ostringstream strstr;
0184     strstr << name_base << " - accumulated thickness = "
0185            << (z_start + geom.get_total_thickness() / 2) / um
0186            << " um expected thickness = " << geom.get_total_thickness() / um
0187            << " um";
0188     G4Exception(
0189         (std::string("PHG4SectorConstructor::Construct_Sectors::") + (name_base)).c_str(),
0190         __FILE__, FatalException, strstr.str().c_str());
0191   }
0192
0193   m_DisplayAction->AddVolume(DetectorLog_Det, "DetectorBox");
0194   if (Verbosity() > 1)
0195   {
0196     std::cout << "PHG4SectorConstructor::Construct_Sectors::" << name_base
0197               << " - total thickness = " << geom.get_total_thickness() / cm << " cm"
0198               << std::endl;
0199     std::cout << "PHG4SectorConstructor::Construct_Sectors::" << name_base << " - "
0200               << map_log_vol.size() << " logical volume constructed" << std::endl;
0201     std::cout << "PHG4SectorConstructor::Construct_Sectors::" << name_base << " - "
0202               << map_phy_vol.size() << " physical volume constructed; "
0203               << map_active_phy_vol.size() << " is active." << std::endl;
0204   }
0205 }
0206
0207 G4VSolid *
0208 PHG4Sector::PHG4SectorConstructor::Construct_Sectors_Plane(  //
0209     const std::string &name,                                 //
0210     const double start_z,                                    //
0211     const double thickness,                                  //
0212     G4VSolid *SecConeBoundary_Det                            //
0213 )
0214 {
0215   assert(SecConeBoundary_Det);
0216
0217   G4VSolid *Sol_Raw = new G4Tubs(name + "_Raw",     // const G4String& pName,
0218                                  0,                 //      G4double pRMin,
0219                                  geom.get_max_R(),  //      G4double pRMax,
0220                                  thickness / 2,     //      G4double pDz,
0221                                  0,                 //      G4double pSPhi,
0222                                  2 * pi             //      G4double pDPhi
0223   );
0224
0225   G4VSolid *Sol_Place = new G4DisplacedSolid(name + "_Place", Sol_Raw, nullptr,
0226                                              G4ThreeVector(0, 0, start_z + thickness / 2));
0227
0228   G4VSolid *Sol = new G4IntersectionSolid(name, Sol_Place,
0229                                           SecConeBoundary_Det);
0230
0231   return Sol;
0232   //  return Sol_Place;
0233 }
0234
0235 void PHG4Sector::Sector_Geometry::SetDefault()
0236 {
0237   N_Sector = 8;
0238   material = "G4_AIR";
0239
0240   //! polar angle for the normal vector
0241   normal_polar_angle = 0;
0242
0243   //! polar angle for edges
0244   min_polar_angle = .1;
0245
0246   //! polar angle for edges
0247   max_polar_angle = .3;
0248
0249   //! distance that detector starts from the normal direction
0250   normal_start = 305 * cm;
0251
0252   min_polar_edge = kConeEdge;
0253
0254   max_polar_edge = kConeEdge;
0255 }
0256
0257 G4LogicalVolume *
0258 PHG4Sector::PHG4SectorConstructor::RegisterLogicalVolume(G4LogicalVolume *v)
0259 {
0260   if (!v)
0261   {
0262     std::cout
0263         << "PHG4SectorConstructor::RegisterVolume - Error - invalid volume!"
0264         << std::endl;
0265     return v;
0266   }
0267   if (map_log_vol.find(v->GetName()) != map_log_vol.end())
0268   {
0269     std::cout << "PHG4SectorConstructor::RegisterVolume - Warning - replacing "
0270               << v->GetName() << std::endl;
0271   }
0272
0273   map_log_vol[v->GetName()] = v;
0274
0275   return v;
0276 }
0277
0278 G4PVPlacement *
0279 PHG4Sector::PHG4SectorConstructor::RegisterPhysicalVolume(G4PVPlacement *v,
0280                                                           const bool active)
0281 {
0282   if (!v)
0283   {
0284     std::cout
0285         << "PHG4SectorConstructor::RegisterPhysicalVolume - Error - invalid volume!"
0286         << std::endl;
0287     return v;
0288   }
0289
0290   phy_vol_idx_t id(v->GetName(), v->GetCopyNo());
0291
0292   if (map_phy_vol.find(id) != map_phy_vol.end())
0293   {
0294     std::cout
0295         << "PHG4SectorConstructor::RegisterPhysicalVolume - Warning - replacing "
0296         << v->GetName() << "[" << v->GetCopyNo() << "]" << std::endl;
0297   }
0298
0299   map_phy_vol[id] = v;
0300
0301   if (active)
0302   {
0303     map_active_phy_vol[id] = v;
0304   }
0305
0306   return v;
0307 }
0308
0309 double
0310 PHG4Sector::Sector_Geometry::get_total_thickness() const
0311 {
0312   double sum = 0;
0313   for (const auto &it : layer_list)
0314   {
0315     sum += it.depth;
0316   }
0317   return sum;
0318 }
0319
0320 double
0321 PHG4Sector::Sector_Geometry::get_max_R() const
0322 {
0323   // Geometry check
0324   assert(std::abs(min_polar_angle - normal_polar_angle) < pi / 2);
0325   assert(std::abs(max_polar_angle - normal_polar_angle) < pi / 2);
0326
0327   if (N_Sector <= 2)
0328   {
0329     // Geometry check
0330     if (cos(min_polar_angle + normal_polar_angle) <= 0)
0331     {
0332       std::stringstream strstr;
0333       strstr << "Geometry check failed. " << std::endl;
0334       strstr << "normal_polar_angle = " << normal_polar_angle << std::endl;
0335       strstr << "min_polar_angle = " << min_polar_angle << std::endl;
0336       strstr << "max_polar_angle = " << max_polar_angle << std::endl;
0337       strstr << "cos(min_polar_angle + normal_polar_angle) = "
0338              << cos(min_polar_angle + normal_polar_angle) << std::endl;
0339
0340       G4Exception("Sector_Geometry::get_max_R", __FILE__, FatalException,
0341                   strstr.str().c_str());
0342     }
0343     if (cos(max_polar_angle + normal_polar_angle) <= 0)
0344     {
0345       std::stringstream strstr;
0346       strstr << "Geometry check failed. " << std::endl;
0347       strstr << "normal_polar_angle = " << normal_polar_angle << std::endl;
0348       strstr << "min_polar_angle = " << min_polar_angle << std::endl;
0349       strstr << "max_polar_angle = " << max_polar_angle << std::endl;
0350       strstr << "cos(max_polar_angle + normal_polar_angle) = "
0351              << cos(max_polar_angle + normal_polar_angle) << std::endl;
0352
0353       G4Exception("Sector_Geometry::get_max_R", __FILE__, FatalException,
0354                   strstr.str().c_str());
0355     }
0356
0357     const double max_tan_angle = std::max(
0358         std::abs(tan(min_polar_angle + normal_polar_angle)),
0359         std::abs(tan(max_polar_angle + normal_polar_angle)));
0360
0361     return (get_normal_start() + get_total_thickness()) * sqrt(1 + max_tan_angle * max_tan_angle);
0362   }
0363   else
0364   {
0365     const double max_angle = std::max(
0366         std::abs(min_polar_angle - normal_polar_angle),
0367         std::abs(max_polar_angle - normal_polar_angle));
0368
0369     return (get_normal_start() + get_total_thickness()) * sqrt(1 + pow(tan(max_angle), 2) + pow(tan(2 * pi / N_Sector), 2)) * 2;
0370   }
0371 }
0372
0373 //! add Entrace window and drift volume
0374 //! Ref: P. Abbon et al. The COMPASS experiment at CERN. Nucl. Instrum. Meth., A577:455
0375 //! 518, 2007. arXiv:hep-ex/0703049, doi:10.1016/j.nima.2007.03.026. 3
0376 void PHG4Sector::Sector_Geometry::AddLayers_DriftVol_COMPASS(const double drift_vol_thickness)
0377 {
0378   //  (drift chamber) is enclosed by two Mylarr [52] cathode foils of 25um thickness,
0379   //  coated with about 10um of graphite,
0380
0381   AddLayer("EntranceWindow", "G4_MYLAR", 25 * um, false, 100);
0382   AddLayer("Cathode", "G4_GRAPHITE", 10 * um, false, 100);
0383   AddLayer("DrftVol", material, drift_vol_thickness, true, 100);
0384 }
0385
0386 //! add HBD GEM to the layer list,
0387 //! Ref: W. Anderson et al. Design, Construction, Operation and Performance of a Hadron
0388 //! Blind Detector for the PHENIX Experiment. Nucl. Instrum. Meth., A646:35 58, 2011.
0389 //! arXiv:1103.4277, doi:10.1016/j.nima.2011.04.015. 3.5.1
0390 void PHG4Sector::Sector_Geometry::AddLayers_HBD_GEM(const int n_GEM_layers)
0391 {
0392   // Internal HBD structure
0393   // From doi:10.1016/j.nima.2011.04.015
0394   // Component Material X0 (cm) Thickness (cm) Area (%) Rad. Length (%)
0395   //  Mesh SS 1.67 0.003 11.5 0.021  <- not used for GEMs trackers
0396   //  AddLayer("Mesh", "Steel",
0397   //          0.003 * cm, false, 11.5);
0398
0399   //  //  GEM frames FR4 17.1 0.15x4 6.5 0.228 <- not used for GEMs trackers
0400   //  AddLayer("Frame0", "G10",
0401   //          0.15 * cm, false, 6.5);
0402
0403   for (int gem = 1; gem <= n_GEM_layers; gem++)
0404   {
0405     std::stringstream sid;
0406     sid << gem;
0407
0408     //  GEM Copper 1.43 0.0005x6 64 0.134
0409     AddLayer(G4String("GEMFrontCu") + G4String(sid.str()), "G4_Cu",
0410              0.0005 * cm, false, 64);
0411
0412     //  GEM Kapton 28.6 0.005x3 64 0.034
0413     AddLayer(G4String("GEMKapton") + G4String(sid.str()), "G4_KAPTON",
0414              0.005 * cm, false, 64);
0415
0416     //  GEM Copper 1.43 0.0005x6 64 0.134
0417     AddLayer(G4String("GEMBackCu") + G4String(sid.str()), "G4_Cu",
0418              0.0005 * cm, false, 64);
0419
0420     //  GEM frames FR4 17.1 0.15x4 6.5 0.228
0421     AddLayer(G4String("Frame") + G4String(sid.str()), "G10", 0.15 * cm, false,
0422              6.5);
0423   }
0424
0425   //  PCB Kapton 28.6 0.005 100 0.017
0426   AddLayer(G4String("PCBKapton"), "G4_KAPTON", 0.005 * cm, false, 100);
0427
0428   //  PCB Copper 1.43 0.0005 80 0.028
0429   AddLayer(G4String("PCBCu"), "G4_Cu", 0.0005 * cm, false, 80);
0430
0431   //  Facesheet FR4 17.1 0.025x2 100 0.292
0432   AddLayer("Facesheet", "G10", 0.025 * 2 * cm, false, 100);
0433
0434   //  Panel core Honeycomb 8170 1.905 100 0.023 <- very thin-X0 stuff, ignore
0435
0436   //  Total vessel 0.82
0437   //  Readout
0438 }
0439
0440 //! add HBD readout,
0441 //! Ref: W. Anderson et al. Design, Construction, Operation and Performance of a Hadron
0442 //! Blind Detector for the PHENIX Experiment. Nucl. Instrum. Meth., A646:35 58, 2011.
0443 //! arXiv:1103.4277, doi:10.1016/j.nima.2011.04.015. 3.5.1
0444 void PHG4Sector::Sector_Geometry::AddLayers_HBD_Readout()
0445 {
0446   //  Readout
0447   //  Readout board FR4/copper 17.1/1.43 0.05/0.001 100 0.367
0448   AddLayer(G4String("ReadoutFR4"), "G10", 0.05 * cm, false, 100);
0449   AddLayer(G4String("ReadoutCu"), "G4_Cu", 0.001 * cm, false, 100);
0450
0451   //  Preamps + sockets Copper 1.43 0.0005 100 0.66
0452   //  Total readout 1.03
0453   AddLayer(G4String("SocketsCu"), "G4_Cu", 0.0005 * cm, false, 100);
0454 }
0455
0456 //! Rough AeroGel detector
0457 //! Ref: T. Iijima et al. A Novel type of proximity focusing RICH counter with multiple
0458 //! refractive index aerogel radiator. Nucl. Instrum. Meth., A548:383 390, 2005. arXiv:
0459 //! physics/0504220, doi:10.1016/j.nima.2005.05.030
0460 void PHG4Sector::Sector_Geometry::AddLayers_AeroGel_ePHENIX(const double radiator_length,
0461                                                             const double expansion_length, std::string radiator)
0462 {
0463   if (radiator == "Default")
0464   {
0465     radiator = "ePHENIX_AeroGel";
0466   }
0467
0468   //  Readout board FR4/copper 17.1/1.43 0.05/0.001 100 0.367
0469   AddLayer("EntranceWindow", "G10", 0.05 * cm, false, 100);
0470   AddLayer("AeroGel", radiator, radiator_length, false, 100);
0471   AddLayer("ExpansionVol", "G4_AIR", expansion_length, false, 100);
0472
0473   // Some readout
0474   AddLayer(G4String("ReadoutFR4"), "G10", 0.05 * cm, false, 100);
0475   AddLayer(G4String("ReadoutCu"), "G4_Cu", 0.001 * cm, false, 100);
0476   AddLayer(G4String("SocketsCu"), "G4_Cu", 0.0005 * cm, false, 100);
0477 }