g4simulation/g4tpc/PHG4TpcDetector.cc

0001 #include "PHG4TpcDetector.h"
0002 #include "PHG4TpcDefs.h"
0003 #include "PHG4TpcDisplayAction.h"
0004
0005 #include <g4detectors/PHG4CellDefs.h>
0006 #include <g4detectors/PHG4TpcGeomv1.h>
0007 #include <g4detectors/PHG4TpcGeomContainer.h>
0008
0009 #include <phparameter/PHParameters.h>
0010
0011 #include <g4main/PHG4Detector.h>       // for PHG4Detector
0012 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0013 #include <g4main/PHG4Subsystem.h>
0014
0015 #include <cdbobjects/CDBTTree.h>
0016
0017 #include <ffamodules/CDBInterface.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/getClass.h>
0023 #include <phool/recoConsts.h>
0024 #include <phool/sphenix_constants.h>
0025
0026 #include <TSystem.h>
0027
0028 #include <Geant4/G4LogicalVolume.hh>
0029 #include <Geant4/G4Material.hh>
0030 #include <Geant4/G4NistManager.hh>
0031 #include <Geant4/G4PVPlacement.hh>
0032 #include <Geant4/G4String.hh>  // for G4String
0033 #include <Geant4/G4SystemOfUnits.hh>
0034 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0035 #include <Geant4/G4Tubs.hh>
0036 #include <Geant4/G4UserLimits.hh>
0037 #include <Geant4/G4VPhysicalVolume.hh>  // for G4VPhysicalVolume
0038
0039 #include <algorithm>  // for max, copy
0040 #include <cassert>
0041 #include <cmath>
0042 #include <cstdlib>   // for exit
0043 #include <iostream>  // for basic_ostream::operator<<
0044 #include <map>       // for map
0045 #include <sstream>
0046
0047 class G4VSolid;
0048 class PHCompositeNode;
0049
0050 //_______________________________________________________________
0051 PHG4TpcDetector::PHG4TpcDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam)
0052   : PHG4Detector(subsys, Node, dnam)
0053   , m_DisplayAction(dynamic_cast<PHG4TpcDisplayAction *>(subsys->GetDisplayAction()))
0054   , m_Params(parameters)
0055   , m_ActiveFlag(m_Params->get_int_param("active"))
0056   , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
0057   , m_InnerCageRadius(m_Params->get_double_param("gas_inner_radius") * cm - m_Params->get_double_param("cage_layer_9_thickness") * cm - m_Params->get_double_param("cage_layer_8_thickness") * cm - m_Params->get_double_param("cage_layer_7_thickness") * cm - m_Params->get_double_param("cage_layer_6_thickness") * cm - m_Params->get_double_param("cage_layer_5_thickness") * cm - m_Params->get_double_param("cage_layer_4_thickness") * cm - m_Params->get_double_param("cage_layer_3_thickness") * cm - m_Params->get_double_param("cage_layer_2_thickness") * cm - m_Params->get_double_param("cage_layer_1_thickness") * cm)
0058   , m_OuterCageRadius(m_Params->get_double_param("gas_outer_radius") * cm + m_Params->get_double_param("cage_layer_9_thickness") * cm + m_Params->get_double_param("cage_layer_8_thickness") * cm + m_Params->get_double_param("cage_layer_7_thickness") * cm + m_Params->get_double_param("cage_layer_6_thickness") * cm + m_Params->get_double_param("cage_layer_5_thickness") * cm + m_Params->get_double_param("cage_layer_4_thickness") * cm + m_Params->get_double_param("cage_layer_3_thickness") * cm + m_Params->get_double_param("cage_layer_2_thickness") * cm + m_Params->get_double_param("cage_layer_1_thickness") * cm)
0059 {
0060 }
0061 //_______________________________________________________________
0062 PHG4TpcDetector::~PHG4TpcDetector()
0063 {
0064   delete m_cdbttree;
0065   delete m_G4UserLimits;
0066 }
0067 //_______________________________________________________________
0068 int PHG4TpcDetector::IsInTpc(G4VPhysicalVolume *volume) const
0069 {
0070   if (m_ActiveFlag)
0071   {
0072     if (m_ActiveVolumeSet.contains(volume))
0073     {
0074       return 1;
0075     }
0076   }
0077   if (m_AbsorberActiveFlag)
0078   {
0079     if (m_AbsorberVolumeSet.contains(volume))
0080     {
0081       return -1;
0082     }
0083   }
0084   return 0;
0085 }
0086
0087 //_______________________________________________________________
0088 void PHG4TpcDetector::ConstructMe(G4LogicalVolume *logicWorld)
0089 {
0090   // create Tpc envelope
0091   // tpc consists of (from inside to gas volume, outside is reversed up to now)
0092   // 1st layer cu
0093   // 2nd layer FR4
0094   // 3rd layer HoneyComb
0095   // 4th layer cu
0096   // 5th layer FR4
0097   // 6th layer Kapton
0098   // 7th layer cu
0099   // 8th layer Kapton
0100   // 9th layer cu
0101
0102   double steplimits = m_Params->get_double_param("steplimits") * cm;
0103   if (std::isfinite(steplimits) && steplimits > 0)
0104   {
0105     m_G4UserLimits = new G4UserLimits(steplimits);
0106   }
0107
0108   G4VSolid *tpc_envelope = new G4Tubs("tpc_envelope", m_InnerCageRadius, m_OuterCageRadius, m_Params->get_double_param("tpc_length") * cm / 2., 0., 2 * M_PI);
0109
0110   recoConsts *rc = recoConsts::instance();
0111   G4LogicalVolume *tpc_envelope_logic = new G4LogicalVolume(tpc_envelope,
0112                                                             GetDetectorMaterial(rc->get_StringFlag("WorldMaterial")),
0113                                                             "tpc_envelope");
0114   m_DisplayAction->AddVolume(tpc_envelope_logic, "TpcEnvelope");
0115
0116   ConstructTpcExternalSupports(logicWorld);
0117
0118   ConstructTpcCageVolume(tpc_envelope_logic);
0119   ConstructTpcGasVolume(tpc_envelope_logic);
0120
0121   new G4PVPlacement(nullptr, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm),
0122                     tpc_envelope_logic, "tpc_envelope",
0123                     logicWorld, false, false, OverlapCheck());
0124
0125   // geometry node
0126   add_geometry_node();
0127 }
0128
0129 void PHG4TpcDetector::CreateTpcGasMixture()
0130 {
0131   G4double density;
0132   G4int ncomponents;
0133   G4int natoms;
0134
0135   G4double tpcGasTemperature = (273.15 + m_Params->get_double_param("TPC_gas_temperature")) * kelvin;
0136   G4double tpcGasPressure = m_Params->get_double_param("TPC_gas_pressure") * atmosphere;
0137
0138   G4Material *CF4 = new G4Material("CF4", density = sphenix_constants::CF4_density * mg / cm3, ncomponents = 2, kStateGas, tpcGasTemperature, tpcGasPressure);
0139   CF4->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), natoms = 1);
0140   CF4->AddElement(G4NistManager::Instance()->FindOrBuildElement("F"), natoms = 4);
0141
0142   G4Material *N2 = new G4Material("N2", density = 1.165 * mg / cm3, ncomponents = 1, kStateGas, tpcGasTemperature, tpcGasPressure);
0143   N2->AddElement(G4NistManager::Instance()->FindOrBuildElement("N"), natoms = 2);
0144
0145   // Create isobutane as only butane is in the standard G4Material list (they have different densities)
0146   // https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Appendix/materialNames.html
0147   G4Material *isobutane = new G4Material("isobutane", density = 2.59 * mg / cm3, ncomponents = 2, kStateGas, tpcGasTemperature, tpcGasPressure);
0148   isobutane->AddElement(G4NistManager::Instance()->FindOrBuildElement("C"), natoms = 4);  // Could use AddElement(PHG4Detector::GetDetectorElement("C", true), 4); instead?
0149   isobutane->AddElement(G4NistManager::Instance()->FindOrBuildElement("H"), natoms = 10);
0150
0151   double Ne_frac = m_Params->get_double_param("Ne_frac");
0152   double Ar_frac = m_Params->get_double_param("Ar_frac");
0153   double CF4_frac = m_Params->get_double_param("CF4_frac");
0154   double N2_frac = m_Params->get_double_param("N2_frac");
0155   double isobutane_frac = m_Params->get_double_param("isobutane_frac");
0156
0157   const double Ne_den = G4NistManager::Instance()->FindOrBuildMaterial("G4_Ne")->GetDensity();
0158   const double Ar_den = G4NistManager::Instance()->FindOrBuildMaterial("G4_Ar")->GetDensity();
0159   const double CF4_den = CF4->GetDensity();
0160   const double N2_den = N2->GetDensity();
0161   const double isobutane_den = isobutane->GetDensity();
0162
0163   const double sphenix_tpc_gas_den = (Ne_den * Ne_frac) + (Ar_den * Ar_frac) + (CF4_den * CF4_frac) + (N2_den * N2_frac) + (isobutane_den * isobutane_frac);
0164
0165   G4Material *sPHENIX_tpc_gas = new G4Material("sPHENIX_TPC_Gas", sphenix_tpc_gas_den, ncomponents = 5, kStateGas);  //, tpcGasTemperature, tpcGasPressure);
0166   sPHENIX_tpc_gas->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_Ne"), (Ne_den * Ne_frac) / sphenix_tpc_gas_den);
0167   sPHENIX_tpc_gas->AddMaterial(G4NistManager::Instance()->FindOrBuildMaterial("G4_Ar"), (Ar_den * Ar_frac) / sphenix_tpc_gas_den);
0168   sPHENIX_tpc_gas->AddMaterial(CF4, (CF4_den * CF4_frac) / sphenix_tpc_gas_den);
0169   sPHENIX_tpc_gas->AddMaterial(N2, (N2_den * N2_frac) / sphenix_tpc_gas_den);
0170   sPHENIX_tpc_gas->AddMaterial(isobutane, (isobutane_den * isobutane_frac) / sphenix_tpc_gas_den);
0171 }
0172
0173 int PHG4TpcDetector::ConstructTpcGasVolume(G4LogicalVolume *tpc_envelope)
0174 {
0175   static std::map<int, std::string> tpcgasvolname =
0176       {{PHG4TpcDefs::North, "tpc_gas_north"},
0177        {PHG4TpcDefs::South, "tpc_gas_south"}};
0178
0179   CreateTpcGasMixture();
0180
0181   // Window / central membrane
0182   double tpc_window_thickness = m_Params->get_double_param("window_thickness") * cm;
0183   double tpc_half_length = (m_Params->get_double_param("tpc_length") * cm - tpc_window_thickness) / 2.;
0184
0185   //'window' (modernly called central membrane only) material is ENIG, not Copper:
0186   // thickness in this recipe are just a ratio.  we set the usual thickness below.
0187   std::vector<double> thickness;
0188   std::vector<std::string> material;
0189   material.emplace_back("G4_Ni");
0190   thickness.push_back(.240 * cm);
0191   material.emplace_back("G4_Au");
0192   thickness.push_back(.008 * cm);
0193   G4Material *temp = nullptr;
0194   temp = GetDetectorMaterial("ENIG", false);
0195   if (temp == nullptr)
0196   {
0197     CreateCompositeMaterial("ENIG", material, thickness);  // see new function below
0198   }
0199
0200   G4VSolid *tpc_window = new G4Tubs("tpc_window", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm, tpc_window_thickness / 2., 0., 2 * M_PI);
0201   // we build our CM surface:
0202   G4LogicalVolume *tpc_window_logic = new G4LogicalVolume(tpc_window,
0203                                                           GetDetectorMaterial("ENIG"),
0204                                                           "tpc_window");
0205   // previously:                                                            GetDetectorMaterial(m_Params->get_string_param("window_surface1_material")),
0206
0207   //  G4VisAttributes *visatt = new G4VisAttributes();
0208   //  visatt->SetVisibility(true);
0209   //  visatt->SetForceSolid(true);
0210   //  visatt->SetColor(PHG4TPCColorDefs::tpc_cu_color);
0211   //  tpc_window_logic->SetVisAttributes(visatt);
0212
0213   m_DisplayAction->AddVolume(tpc_window_logic, "TpcWindow");
0214   G4VPhysicalVolume *tpc_window_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0215                                                          tpc_window_logic, "tpc_window",
0216                                                          tpc_envelope, false, PHG4TpcDefs::Window, OverlapCheck());
0217
0218   m_AbsorberVolumeSet.insert(tpc_window_phys);
0219
0220   // now build the FR4 layer beneath that:
0221   //  Window / central membrane core
0222   double tpc_window_surface1_thickness = m_Params->get_double_param("window_surface1_thickness") * cm;
0223   double tpc_window_surface2_thickness = m_Params->get_double_param("window_surface2_thickness") * cm;
0224   double tpc_window_surface2_core_thickness = tpc_window_thickness - 2 * tpc_window_surface1_thickness;
0225   double tpc_window_core_thickness = tpc_window_surface2_core_thickness - 2 * (tpc_window_surface2_thickness);
0226
0227   G4VSolid *tpc_window_surface2_core =
0228       new G4Tubs("tpc_window_surface2_core", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm,
0229                  tpc_window_surface2_core_thickness / 2., 0., 2 * M_PI);
0230   G4LogicalVolume *tpc_window_surface2_core_logic = new G4LogicalVolume(tpc_window_surface2_core,
0231                                                                         GetDetectorMaterial(m_Params->get_string_param("window_surface2_material")),
0232                                                                         "tpc_window_surface2_core");
0233
0234   m_DisplayAction->AddVolume(tpc_window_surface2_core_logic, "TpcWindow");
0235   //  visatt = new G4VisAttributes();
0236   //  visatt->SetVisibility(true);
0237   //  visatt->SetForceSolid(true);
0238   //  visatt->SetColor(PHG4TPCColorDefs::tpc_pcb_color);
0239   //  tpc_window_surface2_core_logic->SetVisAttributes(visatt);
0240   G4VPhysicalVolume *tpc_window_surface2_core_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0241                                                                        tpc_window_surface2_core_logic, "tpc_window_surface2_core",
0242                                                                        tpc_window_logic, false, PHG4TpcDefs::WindowCore, OverlapCheck());
0243
0244   m_AbsorberVolumeSet.insert(tpc_window_surface2_core_phys);
0245
0246   // and now the honeycomb core:
0247   G4VSolid *tpc_window_core =
0248       new G4Tubs("tpc_window", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm,
0249                  tpc_window_core_thickness / 2., 0., 2 * M_PI);
0250   G4LogicalVolume *tpc_window_core_logic = new G4LogicalVolume(tpc_window_core,
0251                                                                GetDetectorMaterial(m_Params->get_string_param("window_core_material")),
0252                                                                "tpc_window_core");
0253
0254   m_DisplayAction->AddVolume(tpc_window_core_logic, "TpcHoneyComb");
0255   G4VPhysicalVolume *tpc_window_core_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0256                                                               tpc_window_core_logic, "tpc_window_core",
0257                                                               tpc_window_surface2_core_logic, false, PHG4TpcDefs::WindowCore, OverlapCheck());
0258
0259   m_AbsorberVolumeSet.insert(tpc_window_core_phys);
0260
0261   // Gas
0262   G4VSolid *tpc_gas = new G4Tubs("tpc_gas", m_Params->get_double_param("gas_inner_radius") * cm, m_Params->get_double_param("gas_outer_radius") * cm, tpc_half_length / 2., 0., 2 * M_PI);
0263
0264   G4LogicalVolume *tpc_gas_logic = new G4LogicalVolume(tpc_gas,
0265                                                        GetDetectorMaterial("sPHENIX_TPC_Gas"),
0266                                                        "tpc_gas");
0267
0268   tpc_gas_logic->SetUserLimits(m_G4UserLimits);
0269   m_DisplayAction->AddVolume(tpc_gas_logic, "TpcGas");
0270   G4VPhysicalVolume *tpc_gas_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, (tpc_half_length + tpc_window_thickness) / 2.),
0271                                                       tpc_gas_logic, tpcgasvolname[PHG4TpcDefs::North],
0272                                                       tpc_envelope, false, PHG4TpcDefs::North, OverlapCheck());
0273   //  std::cout << "north copy no: " << tpc_gas_phys->GetCopyNo() << std::endl;
0274
0275   m_ActiveVolumeSet.insert(tpc_gas_phys);
0276   tpc_gas_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, -(tpc_half_length + tpc_window_thickness) / 2.),
0277                                    tpc_gas_logic, tpcgasvolname[PHG4TpcDefs::South],
0278                                    tpc_envelope, false, PHG4TpcDefs::South, OverlapCheck());
0279
0280   //  std::cout << "south copy no: " << tpc_gas_phys->GetCopyNo() << std::endl;
0281   m_ActiveVolumeSet.insert(tpc_gas_phys);
0282
0283 #if G4VERSION_NUMBER >= 1033
0284   const G4RegionStore *theRegionStore = G4RegionStore::GetInstance();
0285   G4Region *tpcregion = theRegionStore->GetRegion("REGION_TPCGAS");
0286   tpc_gas_logic->SetRegion(tpcregion);
0287   tpcregion->AddRootLogicalVolume(tpc_gas_logic);
0288 #endif
0289   return 0;
0290 }
0291
0292 int PHG4TpcDetector::ConstructTpcExternalSupports(G4LogicalVolume *logicWorld)
0293 {
0294   // note that these elements are outside the tpc logical volume!
0295
0296   // Two two-inch diam. 304 Stainless Steel solid 'hanger beams' at 32.05" from beam center
0297   // at +/- 41.39 degrees left and right of vertical
0298   // stainless steel: 0.695 iron, 0.190 chromium, 0.095 nickel, 0.020 manganese.
0299   // if we're being pedantic, that is.  But store-bought stainless is probably okay.
0300   // G4Material *StainlessSteel = new G4Material("StainlessSteel",   density = 8.02*g/cm3, 5);
0301   // StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Si"), 0.01);
0302   // StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Mn"), 0.02);
0303   // StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Cr"), 0.19);
0304   // StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Ni"), 0.10);
0305   // StainlessSteel->AddMaterial(matman->FindOrBuildMaterial("G4_Fe"), 0.68);
0306   G4Material *stainlessSteel = GetDetectorMaterial("G4_STAINLESS-STEEL");
0307   double inch = 2.54 * cm;
0308   // double hangerAngle = 41.39 * M_PI / 180.;
0309   // double hangerRadius = 32.05 * inch;
0310   double hangerX = 21.193 * inch;
0311   double hangerY = 24.246 * inch;
0312   double hangerDiameter = 2. * inch;
0313   double extra_length = 0.941 * inch;
0314
0315   G4VSolid *hangerSupport = new G4Tubs("tpc_hanger_support", 0, hangerDiameter / 2., (6 * inch) + extra_length, 0., 2 * M_PI);
0316   G4LogicalVolume *hangerSupportLogic = new G4LogicalVolume(hangerSupport,
0317                                                             stainlessSteel,
0318                                                             "tpc_hanger_support");
0319
0320   double tpc_steel_location = (90 * inch) / 2 - 3 * inch - extra_length / 2.0;
0321
0322   m_DisplayAction->AddVolume(hangerSupportLogic, "TpcHangerSupport");
0323   G4VPhysicalVolume *tpc_hanger_support_phys[4] = {nullptr, nullptr, nullptr, nullptr};
0324   tpc_hanger_support_phys[0] = new G4PVPlacement(nullptr, G4ThreeVector(hangerX, hangerY, tpc_steel_location),
0325                                                  hangerSupportLogic, "tpc_hanger_support_northleft",
0326                                                  logicWorld, false, 0, OverlapCheck());
0327   tpc_hanger_support_phys[1] = new G4PVPlacement(nullptr, G4ThreeVector(-hangerX, hangerY, tpc_steel_location),
0328                                                  hangerSupportLogic, "tpc_hanger_support_northright",
0329                                                  logicWorld, false, 1, OverlapCheck());
0330   tpc_hanger_support_phys[2] = new G4PVPlacement(nullptr, G4ThreeVector(hangerX, hangerY, -tpc_steel_location),
0331                                                  hangerSupportLogic, "tpc_hanger_support_southleft",
0332                                                  logicWorld, false, 2, OverlapCheck());
0333   tpc_hanger_support_phys[3] = new G4PVPlacement(nullptr, G4ThreeVector(-hangerX, hangerY, -tpc_steel_location),
0334                                                  hangerSupportLogic, "tpc_hanger_support_southright",
0335                                                  logicWorld, false, 3, OverlapCheck());
0336
0337   m_AbsorberVolumeSet.insert(tpc_hanger_support_phys[0]);
0338   m_AbsorberVolumeSet.insert(tpc_hanger_support_phys[1]);
0339   m_AbsorberVolumeSet.insert(tpc_hanger_support_phys[2]);
0340   m_AbsorberVolumeSet.insert(tpc_hanger_support_phys[3]);
0341
0342   G4Material *carbonFiber = GetDetectorMaterial("CFRP_INTT");
0343   double hangerBeamInnerDiameter = 2.0 * inch;
0344   double hangerBeamOuterDiameter = 2.521 * inch;
0345
0346   G4VSolid *hangerBeam = new G4Tubs("tpc_hanger_beam", hangerBeamInnerDiameter / 2., hangerBeamOuterDiameter / 2., (90 * inch) / 2., 0., 2 * M_PI);
0347   G4LogicalVolume *hangerBeamLogic = new G4LogicalVolume(hangerBeam,
0348                                                          carbonFiber,
0349                                                          "tpc_hanger_beam");
0350
0351   m_DisplayAction->AddVolume(hangerBeamLogic, "TpcHangerBeam");
0352   G4VPhysicalVolume *tpc_hanger_beam_phys[2] = {nullptr, nullptr};
0353   tpc_hanger_beam_phys[0] = new G4PVPlacement(nullptr, G4ThreeVector(hangerX, hangerY, 0),
0354                                               hangerBeamLogic, "tpc_hanger_beam_left",
0355                                               logicWorld, false, 0, OverlapCheck());
0356   tpc_hanger_beam_phys[1] = new G4PVPlacement(nullptr, G4ThreeVector(-hangerX, hangerY, 0),
0357                                               hangerBeamLogic, "tpc_hanger_beam_right",
0358                                               logicWorld, false, 1, OverlapCheck());
0359
0360   m_AbsorberVolumeSet.insert(tpc_hanger_beam_phys[0]);
0361   m_AbsorberVolumeSet.insert(tpc_hanger_beam_phys[1]);
0362
0363   // Twelve one-inch diam carbon fiber rods of thickness 1/16" at 31.04" from beam center
0364   // borrowed from the INTT specification of carbon fiber
0365   // note that this defines a clocking!
0366   double rodAngleStart = M_PI / 12.;
0367   double rodAngularSpacing = 2 * M_PI / 12.;
0368   double rodRadius = 31.5 * inch;
0369   double rodWallThickness = 1. / 8. * inch;
0370   double rodDiameter = 3. / 4. * inch;
0371   G4VSolid *tieRod = new G4Tubs("tpc_tie_rod", rodDiameter / 2. - rodWallThickness, rodDiameter / 2., (m_Params->get_double_param("tpc_length") * cm) / 2., 0., 2 * M_PI);
0372   G4LogicalVolume *tieRodLogic = new G4LogicalVolume(tieRod,
0373                                                      carbonFiber,
0374                                                      "tpc_tie_rod");
0375
0376   G4VPhysicalVolume *tpc_tie_rod_phys[12] = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr,
0377                                              nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
0378
0379   std::ostringstream name;
0380   for (int i = 0; i < 12; i++)
0381   {
0382     double ang = rodAngleStart + rodAngularSpacing * i;
0383     name.str("");
0384     name << "tpc_tie_rod_" << i;
0385     tpc_tie_rod_phys[i] = new G4PVPlacement(nullptr, G4ThreeVector(rodRadius * cos(ang), rodRadius * sin(ang), 0),
0386                                             tieRodLogic, name.str(),
0387                                             logicWorld, false, i, OverlapCheck());
0388     m_AbsorberVolumeSet.insert(tpc_tie_rod_phys[i]);
0389   }
0390
0391   //  G4VisAttributes *visatt = new G4VisAttributes();
0392   //  visatt->SetVisibility(true);
0393   //  visatt->SetForceSolid(true);
0394   //  visatt->SetColor(PHG4TPCColorDefs::tpc_cu_color);
0395   //  tpc_window_logic->SetVisAttributes(visatt);
0396
0397   return 0;
0398 }
0399
0400 int PHG4TpcDetector::ConstructTpcCageVolume(G4LogicalVolume *tpc_envelope)
0401 {
0402   // 8th layer cu
0403   // 9th layer FR4
0404   // 10th layer HoneyComb
0405   // 11th layer FR4
0406   // 12th layer Kapton
0407   // 13th layer FR4
0408   static const int nlayers = 9;
0409   static const double thickness[nlayers] = {m_Params->get_double_param("cage_layer_1_thickness") * cm,
0410                                             m_Params->get_double_param("cage_layer_2_thickness") * cm,
0411                                             m_Params->get_double_param("cage_layer_3_thickness") * cm,
0412                                             m_Params->get_double_param("cage_layer_4_thickness") * cm,
0413                                             m_Params->get_double_param("cage_layer_5_thickness") * cm,
0414                                             m_Params->get_double_param("cage_layer_6_thickness") * cm,
0415                                             m_Params->get_double_param("cage_layer_7_thickness") * cm,
0416                                             m_Params->get_double_param("cage_layer_8_thickness") * cm,
0417                                             m_Params->get_double_param("cage_layer_9_thickness") * cm};
0418
0419   static const std::string material[nlayers] = {m_Params->get_string_param("cage_layer_1_material"),
0420                                                 m_Params->get_string_param("cage_layer_2_material"),
0421                                                 m_Params->get_string_param("cage_layer_3_material"),
0422                                                 m_Params->get_string_param("cage_layer_4_material"),
0423                                                 m_Params->get_string_param("cage_layer_5_material"),
0424                                                 m_Params->get_string_param("cage_layer_6_material"),
0425                                                 m_Params->get_string_param("cage_layer_7_material"),
0426                                                 m_Params->get_string_param("cage_layer_8_material"),
0427                                                 m_Params->get_string_param("cage_layer_9_material")};
0428
0429   double tpc_cage_radius = m_InnerCageRadius;
0430   std::ostringstream name;
0431   for (int i = 0; i < nlayers; i++)
0432   {
0433     name.str("");
0434     int layerno = i + 1;
0435     name << "tpc_cage_layer_" << layerno;
0436     G4VSolid *tpc_cage_layer = new G4Tubs(name.str(), tpc_cage_radius, tpc_cage_radius + thickness[i], m_Params->get_double_param("tpc_length") * cm / 2., 0., 2 * M_PI);
0437     G4LogicalVolume *tpc_cage_layer_logic = new G4LogicalVolume(tpc_cage_layer,
0438                                                                 GetDetectorMaterial(material[i]),
0439                                                                 name.str());
0440     m_DisplayAction->AddTpcInnerLayer(tpc_cage_layer_logic);
0441     G4VPhysicalVolume *tpc_cage_layer_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0442                                                                tpc_cage_layer_logic, name.str(),
0443                                                                tpc_envelope, false, layerno, OverlapCheck());
0444     m_AbsorberVolumeSet.insert(tpc_cage_layer_phys);
0445     tpc_cage_radius += thickness[i];
0446   }
0447   // outer cage
0448
0449   tpc_cage_radius = m_OuterCageRadius;
0450   for (int i = 0; i < nlayers; i++)
0451   {
0452     tpc_cage_radius -= thickness[i];
0453     name.str("");
0454     int layerno = 10 + 1 + i;  // so the accompanying inner layer is layer - 10
0455     name << "tpc_cage_layer_" << layerno;
0456     G4VSolid *tpc_cage_layer = new G4Tubs(name.str(), tpc_cage_radius, tpc_cage_radius + thickness[i], m_Params->get_double_param("tpc_length") * cm / 2., 0., 2 * M_PI);
0457     G4LogicalVolume *tpc_cage_layer_logic = new G4LogicalVolume(tpc_cage_layer,
0458                                                                 GetDetectorMaterial(material[i]),
0459                                                                 name.str());
0460     m_DisplayAction->AddTpcOuterLayer(tpc_cage_layer_logic);
0461     G4VPhysicalVolume *tpc_cage_layer_phys = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0),
0462                                                                tpc_cage_layer_logic, name.str(),
0463                                                                tpc_envelope, false, layerno, OverlapCheck());
0464     m_AbsorberVolumeSet.insert(tpc_cage_layer_phys);
0465   }
0466
0467   return 0;
0468 }
0469
0470 void PHG4TpcDetector::CreateCompositeMaterial(
0471     const std::string &compositeName,
0472     std::vector<std::string> materialName,
0473     const std::vector<double> &thickness)
0474 {
0475   // takes in a list of material names known to Geant already, and thicknesses, and creates a new material called compositeName.
0476
0477   // check that desired material name doesn't already exist
0478   G4Material *tempmat = GetDetectorMaterial(compositeName, false);
0479
0480   if (tempmat != nullptr)
0481   {
0482     std::cout << __PRETTY_FUNCTION__ << " Fatal Error: composite material " << compositeName << " already exists" << std::endl;
0483     assert(!tempmat);
0484   }
0485
0486   // check that both arrays have the same depth
0487   assert(materialName.size() == thickness.size());
0488
0489   // sum up the areal density and total thickness so we can divvy it out
0490   double totalArealDensity = 0;
0491   double totalThickness = 0;
0492   for (std::vector<double>::size_type i = 0; i < thickness.size(); i++)
0493   {
0494     tempmat = GetDetectorMaterial(materialName[i]);
0495     if (tempmat == nullptr)
0496     {
0497       std::cout << __PRETTY_FUNCTION__ << " Fatal Error: component material " << materialName[i] << " does not exist." << std::endl;
0498       gSystem->Exit(1);
0499       exit(1);
0500     }
0501     totalArealDensity += tempmat->GetDensity() * thickness[i];
0502     totalThickness += thickness[i];
0503   }
0504
0505   // register a new material with the average density of the whole:
0506   double compositeDensity = totalArealDensity / totalThickness;
0507   G4Material *composite = new G4Material(compositeName, compositeDensity, thickness.size());
0508
0509   // now calculate the fraction due to each material, and register those
0510   for (std::vector<double>::size_type i = 0; i < thickness.size(); i++)
0511   {
0512     tempmat = GetDetectorMaterial(materialName[i]);  // don't need to check this, since we did in the previous loop.
0513     composite->AddMaterial(tempmat, thickness[i] * tempmat->GetDensity() / totalArealDensity);
0514   }
0515
0516   // how to register our finished material?
0517   return;
0518 }
0519
0520 //_______________________________________________________________
0521 void PHG4TpcDetector::add_geometry_node()
0522 {
0523   // create PHG4TpcGeomContainer and put on node tree
0524   const std::string geonode_name = "TPCGEOMCONTAINER";
0525   auto *geonode = findNode::getClass<PHG4TpcGeomContainer>(topNode(), geonode_name);
0526   if (!geonode)
0527   {
0528     geonode = new PHG4TpcGeomContainer;
0529     PHNodeIterator iter(topNode());
0530     auto *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0531     PHNodeIterator runiter(runNode);
0532     PHCompositeNode *geomNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", "RECO_TRACKING_GEOMETRY"));
0533     if(!geomNode)
0534       {
0535         geomNode = new PHCompositeNode("RECO_TRACKING_GEOMETRY");
0536         runNode->addNode(geomNode);
0537       }
0538     auto *newNode = new PHIODataNode<PHObject>(geonode, geonode_name, "PHObject");
0539     geomNode->addNode(newNode);
0540   }
0541
0542   m_cdb = CDBInterface::instance();
0543   std::string calibdir = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP");
0544
0545   if (!calibdir.empty())
0546   {
0547     m_cdbttree = new CDBTTree(calibdir);
0548     m_cdbttree->LoadCalibrations();
0549   }
0550   else
0551   {
0552     std::cout << "PHG4TpcPadPlaneReadout::InitRun no TPC_FEE_CHANNEL_MAP found" << std::endl;
0553     exit(1);
0554   }
0555
0556   const std::array<int, 3> NTpcLayers =
0557       {{m_Params->get_int_param("ntpc_layers_inner"),
0558         m_Params->get_int_param("ntpc_layers_mid"),
0559         m_Params->get_int_param("ntpc_layers_outer")}};
0560
0561   std::array<double, 3> MinLayer{};
0562   MinLayer[0] = m_Params->get_int_param("tpc_minlayer_inner");
0563   MinLayer[1] = MinLayer[0] + NTpcLayers[0];
0564   MinLayer[2] = MinLayer[1] + NTpcLayers[1];
0565
0566   const std::array<double, 5> Thickness =
0567       {{
0568           0.56598621677629212,
0569           1.0206889851687158,
0570           1.0970475085472556,
0571           0.5630547309825637,
0572           0.56891770257002054,
0573       }};
0574
0575   const double drift_velocity_sim = m_Params->get_double_param("drift_velocity_sim");
0576   const double tpc_adc_clock = m_Params->get_double_param("tpc_adc_clock");
0577   const double maxdriftlength = m_Params->get_double_param("maxdriftlength");
0578   const double TBinWidth = tpc_adc_clock;
0579   const double extended_readout_time = m_Params->get_double_param("extended_readout_time");
0580   const double MaxT = extended_readout_time + 2. * maxdriftlength / drift_velocity_sim;  // allows for extended time readout
0581   const double MinT = 0;
0582   const int NTBins = (int) ((MaxT - MinT) / TBinWidth) + 1;
0583
0584   std::cout << PHWHERE << "MaxT " << MaxT << " TBinWidth " << TBinWidth << " extended readout time "
0585             << extended_readout_time << " NTBins = " << NTBins << " drift velocity sim " << drift_velocity_sim << std::endl;
0586
0587   const std::array<int, 3> NPhiBins =
0588       {{m_Params->get_int_param("ntpc_phibins_inner"),
0589         m_Params->get_int_param("ntpc_phibins_mid"),
0590         m_Params->get_int_param("ntpc_phibins_outer")}};
0591
0592   // should move to a common file
0593   static constexpr int NSides = 2;
0594   static constexpr int NSectors = 12;
0595   static constexpr int NLayers = 16 * 3;
0596
0597   std::array<std::vector<double>, NSides> sector_R_bias;
0598   std::array<std::vector<double>, NSides> sector_Phi_bias;
0599   std::array<std::vector<double>, NSides> sector_min_Phi;
0600   std::array<std::vector<double>, NSides> sector_max_Phi;
0601   std::array<std::vector<double>, NLayers> pad_phi;
0602   std::array<std::vector<double>, NLayers> pad_R;
0603   std::array<double, NLayers> layer_radius{};
0604   std::array<double, NLayers> phi_bin_width_cdb{};
0605   std::array<std::array<std::array<double, 3>, NSectors>, NSides> sec_max_phi{};
0606   std::array<std::array<std::array<double, 3>, NSectors>, NSides> sec_min_phi{};
0607   int Nfee = 26;
0608   int Nch = 256;
0609
0610   for (int f = 0; f < Nfee; f++)
0611   {
0612     for (int ch = 0; ch < Nch; ch++)
0613     {
0614       unsigned int key = 256 * (f) + ch;
0615       std::string varname = "layer";
0616       int l = m_cdbttree->GetIntValue(key, varname);
0617       if (l > 6)
0618       {
0619         int v_layer = l - 7;
0620         std::string phiname = "phi";
0621         pad_phi[v_layer].push_back(m_cdbttree->GetDoubleValue(key, phiname));
0622         std::string rname = "R";
0623         pad_R[v_layer].push_back(m_cdbttree->GetDoubleValue(key, rname));
0624       }
0625     }
0626   }
0627
0628   for (size_t layer = 0; layer < NLayers; layer++)
0629   {
0630     layer_radius[layer] = 0;
0631     for (int pad = 0; pad < (int) pad_R[(int) layer].size(); pad++)
0632     {
0633       layer_radius[(int) layer] += pad_R[(int) layer][pad];
0634     }
0635
0636     layer_radius[(int) layer] = layer_radius[(int) layer] / pad_R[(int) layer].size();
0637     layer_radius[(int) layer] = layer_radius[(int) layer] / 10.;
0638
0639     auto min_phi_iter = std::min_element(pad_phi[layer].begin(), pad_phi[layer].end());
0640     auto max_phi_iter = std::max_element(pad_phi[layer].begin(), pad_phi[layer].end());
0641     double min_phi = static_cast<double>(*min_phi_iter);
0642     double max_phi = static_cast<double>(*max_phi_iter);
0643     min_phi = min_phi - M_PI / 2.;
0644     max_phi = max_phi - M_PI / 2.;
0645     phi_bin_width_cdb[layer] = std::abs(max_phi - min_phi) / (NPhiBins[(int) (layer / 16)] / 12 - 1);  // NOLINT(bugprone-integer-division)
0646     double SectorPhi = std::abs(max_phi - min_phi) + phi_bin_width_cdb[layer];
0647     for (int zside = 0; zside < 2; zside++)
0648     {
0649       for (int isector = 0; isector < NSectors; isector++)  // 12 sectors
0650       {
0651         if (zside == 0)
0652         {
0653           sec_min_phi[zside][isector][(int) layer / 16] = M_PI - 2 * M_PI / 12 * (isector + 1) + (-(max_phi) -phi_bin_width_cdb[layer] / 2.);
0654           sec_max_phi[zside][isector][(int) layer / 16] = sec_min_phi[zside][isector][(int) layer / 16] + SectorPhi;
0655         }
0656         if (zside == 1)
0657         {
0658           sec_max_phi[zside][isector][(int) layer / 16] = M_PI - 2 * M_PI / 12 * (isector + 1) + (max_phi + phi_bin_width_cdb[layer] / 2.);
0659           sec_min_phi[zside][isector][(int) layer / 16] = sec_max_phi[zside][isector][(int) layer / 16] - SectorPhi;
0660         }
0661       }
0662     }
0663   }
0664
0665   for (int iregion = 0; iregion < 3; ++iregion)
0666   {
0667     // int zside = 0;
0668     for (int zside = 0; zside < 2; zside++)
0669     {
0670       sector_R_bias[zside].clear();
0671       sector_Phi_bias[zside].clear();
0672       sector_min_Phi[zside].clear();
0673       sector_max_Phi[zside].clear();
0674       // int eff_layer = 0;
0675       for (int isector = 0; isector < NSectors; ++isector)  // 12 sectors
0676       {
0677         // no bias per default
0678         // TODO: confirm with what is in PHG4TpcPadPlane Readout
0679         sector_R_bias[zside].push_back(0);
0680         sector_Phi_bias[zside].push_back(0);
0681
0682         sector_min_Phi[zside].push_back(sec_min_phi[zside][isector][iregion]);
0683         sector_max_Phi[zside].push_back(sec_max_phi[zside][isector][iregion]);
0684       }  // isector
0685     }
0686
0687     for (int layer = MinLayer[iregion]; layer < MinLayer[iregion] + NTpcLayers[iregion]; ++layer)
0688     {
0689       if (Verbosity())
0690       {
0691         std::cout << " layer " << layer << " MinLayer " << MinLayer[iregion] << " region " << iregion
0692                   << " radius " << layer_radius[layer - 7]
0693                   << " thickness " << Thickness[iregion]
0694                   << " NTBins " << NTBins << " tmin " << MinT << " tstep " << TBinWidth
0695                   << " phibins " << NPhiBins[iregion] << " phistep " << phi_bin_width_cdb[layer] << std::endl;
0696       }
0697
0698       auto *layerseggeo = new PHG4TpcGeomv1;
0699       layerseggeo->set_layer(layer);
0700
0701       double r_length = Thickness[iregion];
0702       if (iregion == 0)
0703       {
0704         if (layer % 2 == 0)
0705         {
0706           r_length = Thickness[4];
0707         }
0708         else
0709         {
0710           r_length = Thickness[3];
0711         }
0712       }
0713       int v_layer = layer - 7;
0714       if (v_layer >= 0)
0715       {
0716         layerseggeo->set_thickness(r_length);
0717         layerseggeo->set_radius(layer_radius[v_layer]);
0718         layerseggeo->set_binning(PHG4CellDefs::sizebinning);
0719         layerseggeo->set_zbins(NTBins);
0720         layerseggeo->set_zmin(MinT);
0721         layerseggeo->set_zstep(TBinWidth);
0722         layerseggeo->set_phibins(NPhiBins[iregion]);
0723         layerseggeo->set_phistep(phi_bin_width_cdb[v_layer]);
0724         layerseggeo->set_r_bias(sector_R_bias);
0725         layerseggeo->set_phi_bias(sector_Phi_bias);
0726         layerseggeo->set_sector_min_phi(sector_min_Phi);
0727         layerseggeo->set_sector_max_phi(sector_max_Phi);
0728     // capture the z geometry setup in the layergeom
0729     layerseggeo->set_max_driftlength(m_Params->get_double_param("maxdriftlength"));
0730     layerseggeo->set_CM_halfwidth(m_Params->get_double_param("CM_halfwidth"));
0731     layerseggeo->set_adc_clock(m_Params->get_double_param("tpc_adc_clock"));
0732     layerseggeo->set_extended_readout_time(m_Params->get_double_param("extended_readout_time"));
0733     layerseggeo->set_drift_velocity_sim(m_Params->get_double_param("drift_velocity_sim"));
0734       }
0735
0736       // Chris Pinkenburg: greater causes huge memory growth which causes problems
0737       // on our farm. If you need to increase this - TALK TO ME first
0738       if (NPhiBins[iregion] * NTBins > 5100000)
0739       {
0740         std::cout << "increase Tpc cellsize, number of cells "
0741                   << NPhiBins[iregion] * NTBins << " for layer " << layer
0742                   << " exceed 5.1M limit" << std::endl;
0743         gSystem->Exit(1);
0744       }
0745       geonode->AddLayerCellGeom(layerseggeo);
0746     }
0747   }
0748 }