Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:18:16 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!