File indexing completed on 2025-08-06 08:19:16
0001 #include "PHG4IHCalDetector.h"
0002
0003 #include "PHG4IHCalDisplayAction.h"
0004
0005 #include <g4detectors/PHG4DetectorSubsystem.h>
0006 #include <g4detectors/PHG4HcalDefs.h>
0007
0008 #include <phparameter/PHParameters.h>
0009
0010 #include <g4gdml/PHG4GDMLConfig.hh>
0011 #include <g4gdml/PHG4GDMLUtility.hh>
0012
0013 #include <calobase/RawTowerDefs.h> // for convert_name_...
0014 #include <calobase/RawTowerGeom.h> // for RawTowerGeom
0015 #include <calobase/RawTowerGeomContainer.h> // for RawTowerGeomC...
0016 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0017 #include <calobase/RawTowerGeomv1.h>
0018
0019 #include <g4main/PHG4Detector.h>
0020 #include <g4main/PHG4DisplayAction.h>
0021 #include <g4main/PHG4Subsystem.h>
0022 #include <g4main/PHG4Utils.h>
0023
0024 #include <ffamodules/CDBInterface.h>
0025
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/PHIODataNode.h>
0028 #include <phool/PHNode.h> // for PHNode
0029 #include <phool/PHNodeIterator.h>
0030 #include <phool/PHObject.h> // for PHObject
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>
0033 #include <phool/recoConsts.h>
0034
0035 #include <TSystem.h>
0036
0037 #include <Geant4/G4AssemblyVolume.hh>
0038 #include <Geant4/G4IonisParamMat.hh>
0039 #include <Geant4/G4LogicalVolume.hh>
0040 #include <Geant4/G4Material.hh>
0041 #include <Geant4/G4MaterialTable.hh>
0042 #include <Geant4/G4PVPlacement.hh>
0043 #include <Geant4/G4RotationMatrix.hh>
0044 #include <Geant4/G4String.hh>
0045 #include <Geant4/G4SystemOfUnits.hh>
0046 #include <Geant4/G4ThreeVector.hh>
0047 #include <Geant4/G4Transform3D.hh>
0048 #include <Geant4/G4Tubs.hh>
0049 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
0050 #include <Geant4/G4VSolid.hh>
0051
0052 #pragma GCC diagnostic push
0053 #pragma GCC diagnostic ignored "-Wshadow"
0054 #pragma GCC diagnostic ignored "-Wpedantic"
0055 #include <Geant4/G4GDMLParser.hh>
0056 #include <Geant4/G4GDMLReadStructure.hh> // for G4GDMLReadStructure
0057 #pragma GCC diagnostic pop
0058
0059 #include <boost/lexical_cast.hpp>
0060 #include <boost/tokenizer.hpp>
0061
0062 #include <cassert>
0063 #include <cstdlib>
0064 #include <cmath>
0065 #include <filesystem>
0066 #include <iostream>
0067 #include <memory> // for unique_ptr
0068 #include <utility> // for pair, make_pair
0069 #include <vector> // for vector, vector<>::iter...
0070
0071 PHG4IHCalDetector::PHG4IHCalDetector(PHG4Subsystem *subsys, PHCompositeNode *Node, PHParameters *parameters, const std::string &dnam)
0072 : PHG4Detector(subsys, Node, dnam)
0073 , m_DisplayAction(dynamic_cast<PHG4IHCalDisplayAction *>(subsys->GetDisplayAction()))
0074 , m_Params(parameters)
0075 , m_InnerRadius(m_Params->get_double_param("inner_radius") * cm)
0076 , m_OuterRadius(m_Params->get_double_param("outer_radius") * cm)
0077 , m_SizeZ(m_Params->get_double_param("size_z") * cm)
0078 , m_NumScintiPlates(m_Params->get_int_param(PHG4HcalDefs::scipertwr) * m_Params->get_int_param("n_towers"))
0079 , m_Active(m_Params->get_int_param("active"))
0080 , m_AbsorberActive(m_Params->get_int_param("absorberactive"))
0081 , m_GDMPath(m_Params->get_string_param("GDMPath"))
0082 {
0083 gdml_config = PHG4GDMLUtility::GetOrMakeConfigNode(Node);
0084 assert(gdml_config);
0085
0086
0087 if (std::filesystem::path(m_GDMPath).extension() != ".gdml")
0088 {
0089 m_GDMPath = CDBInterface::instance()->getUrl(m_GDMPath);
0090 m_Params->set_string_param("GDMPath", m_GDMPath);
0091 }
0092 }
0093
0094 PHG4IHCalDetector::~PHG4IHCalDetector()
0095 {
0096 delete m_ScintiMotherAssembly;
0097 }
0098
0099
0100
0101 int PHG4IHCalDetector::IsInIHCal(G4VPhysicalVolume *volume) const
0102 {
0103 if (m_AbsorberActive)
0104 {
0105 if (m_SteelAbsorberLogVolSet.find(volume->GetLogicalVolume()) != m_SteelAbsorberLogVolSet.end())
0106 {
0107 return -1;
0108 }
0109 }
0110 if (m_Active)
0111 {
0112 if (m_ScintiTileLogVolSet.find(volume->GetLogicalVolume()) != m_ScintiTileLogVolSet.end())
0113 {
0114 return 1;
0115 }
0116 }
0117 return 0;
0118 }
0119
0120
0121
0122 void PHG4IHCalDetector::ConstructMe(G4LogicalVolume *logicWorld)
0123 {
0124 recoConsts *rc = recoConsts::instance();
0125 G4Material *worldmat = GetDetectorMaterial(rc->get_StringFlag("WorldMaterial"));
0126 G4VSolid *hcal_envelope_cylinder = new G4Tubs("IHCal_envelope_solid", m_InnerRadius, m_OuterRadius, m_SizeZ / 2., 0, 2 * M_PI);
0127 m_VolumeEnvelope = hcal_envelope_cylinder->GetCubicVolume();
0128 G4LogicalVolume *hcal_envelope_log = new G4LogicalVolume(hcal_envelope_cylinder, worldmat, "Hcal_envelope", nullptr, nullptr, nullptr);
0129
0130 G4RotationMatrix hcal_rotm;
0131 hcal_rotm.rotateX(m_Params->get_double_param("rot_x") * deg);
0132 hcal_rotm.rotateY(m_Params->get_double_param("rot_y") * deg);
0133 hcal_rotm.rotateZ(m_Params->get_double_param("rot_z") * deg);
0134 G4VPhysicalVolume *mothervol = new G4PVPlacement(G4Transform3D(hcal_rotm, 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)), hcal_envelope_log, "IHCalEnvelope", logicWorld, false, false, OverlapCheck());
0135 m_DisplayAction->SetMyTopVolume(mothervol);
0136 ConstructIHCal(hcal_envelope_log);
0137
0138
0139 assert(gdml_config);
0140 gdml_config->exclude_physical_vol(mothervol);
0141 gdml_config->exclude_logical_vol(hcal_envelope_log);
0142
0143 const G4MaterialTable *mtable = G4Material::GetMaterialTable();
0144 int nMaterials = G4Material::GetNumberOfMaterials();
0145 for (auto i = 0; i < nMaterials; ++i)
0146 {
0147 const G4Material *mat = (*mtable)[i];
0148 if (mat->GetName() == "Uniplast_scintillator")
0149 {
0150 if ((mat->GetIonisation()->GetBirksConstant()) == 0)
0151 {
0152 mat->GetIonisation()->SetBirksConstant(m_Params->get_double_param("Birk_const"));
0153 }
0154 }
0155 }
0156 if (!m_Params->get_int_param("saveg4hit"))
0157 {
0158 AddGeometryNode();
0159 }
0160 return;
0161 }
0162
0163 int PHG4IHCalDetector::ConstructIHCal(G4LogicalVolume *hcalenvelope)
0164 {
0165
0166 std::unique_ptr<G4GDMLReadStructure> reader(new G4GDMLReadStructure());
0167 G4GDMLParser gdmlParser(reader.get());
0168 gdmlParser.SetOverlapCheck(OverlapCheck());
0169 if (!std::filesystem::exists(m_GDMPath))
0170 {
0171 std::cout << PHWHERE << " Inner HCal gdml file " << m_GDMPath << " not found" << std::endl;
0172 gSystem->Exit(1);
0173 exit(1);
0174 }
0175 gdmlParser.Read(m_GDMPath, false);
0176
0177 G4AssemblyVolume *abs_asym = reader->GetAssembly("InnerSector");
0178 m_ScintiMotherAssembly = reader->GetAssembly("InnerTileAssembly90");
0179 std::vector<G4VPhysicalVolume *>::iterator it = abs_asym->GetVolumesIterator();
0180 static const unsigned int tilepersec = 24 * 4 * 2;
0181 for (unsigned int isector = 0; isector < abs_asym->TotalImprintedVolumes(); isector++)
0182 {
0183 m_DisplayAction->AddSteelVolume((*it)->GetLogicalVolume());
0184 m_SteelAbsorberLogVolSet.insert((*it)->GetLogicalVolume());
0185 hcalenvelope->AddDaughter((*it));
0186 m_AbsorberPhysVolMap.insert(std::make_pair(*it, isector));
0187 m_VolumeSteel += (*it)->GetLogicalVolume()->GetSolid()->GetCubicVolume();
0188 std::vector<G4VPhysicalVolume *>::iterator its = m_ScintiMotherAssembly->GetVolumesIterator();
0189 unsigned int ioff = isector * tilepersec;
0190 for (unsigned int j = 0; j < ioff; j++)
0191 {
0192 ++its;
0193 }
0194 for (unsigned int j = ioff; j < ioff + tilepersec; j++)
0195 {
0196 m_DisplayAction->AddScintiVolume((*its)->GetLogicalVolume());
0197 m_ScintiTileLogVolSet.insert((*its)->GetLogicalVolume());
0198 hcalenvelope->AddDaughter((*its));
0199 m_ScintiTilePhysVolMap.insert(std::make_pair(*its, ExtractLayerTowerId(isector, *its)));
0200 m_VolumeScintillator += (*its)->GetLogicalVolume()->GetSolid()->GetCubicVolume();
0201 ++its;
0202 }
0203
0204 ++it;
0205 }
0206 return 0;
0207 }
0208
0209 int PHG4IHCalDetector::ConsistencyCheck() const
0210 {
0211
0212 if (m_InnerRadius >= m_OuterRadius)
0213 {
0214 std::cout << PHWHERE << ": Inner Radius " << m_InnerRadius / cm
0215 << " cm larger than Outer Radius " << m_OuterRadius / cm
0216 << " cm" << std::endl;
0217 gSystem->Exit(1);
0218 }
0219 return 0;
0220 }
0221
0222 void PHG4IHCalDetector::Print(const std::string &what) const
0223 {
0224 std::cout << "Inner Hcal Detector:" << std::endl;
0225 if (what == "ALL" || what == "VOLUME")
0226 {
0227 std::cout << "Volume Envelope: " << m_VolumeEnvelope / cm3 << " cm^3" << std::endl;
0228 std::cout << "Volume Steel: " << m_VolumeSteel / cm3 << " cm^3" << std::endl;
0229 std::cout << "Volume Scintillator: " << m_VolumeScintillator / cm3 << " cm^3" << std::endl;
0230 std::cout << "Volume Air: " << (m_VolumeEnvelope - m_VolumeSteel - m_VolumeScintillator) / cm3 << " cm^3" << std::endl;
0231 }
0232 return;
0233 }
0234
0235 std::tuple<int, int, int> PHG4IHCalDetector::GetLayerTowerId(G4VPhysicalVolume *volume) const
0236 {
0237 auto it = m_ScintiTilePhysVolMap.find(volume);
0238 if (it != m_ScintiTilePhysVolMap.end())
0239 {
0240 return it->second;
0241 }
0242 std::cout << "could not locate volume " << volume->GetName()
0243 << " in Inner Hcal scintillator map" << std::endl;
0244 gSystem->Exit(1);
0245
0246
0247 exit(1);
0248 }
0249
0250 int PHG4IHCalDetector::GetSectorId(G4VPhysicalVolume *volume) const
0251 {
0252 auto it = m_AbsorberPhysVolMap.find(volume);
0253 if (it != m_AbsorberPhysVolMap.end())
0254 {
0255 return it->second;
0256 }
0257 std::cout << "could not locate volume " << volume->GetName()
0258 << " in Inner Hcal Absorber map" << std::endl;
0259 gSystem->Exit(1);
0260
0261
0262 exit(1);
0263 }
0264
0265 std::tuple<int, int, int> PHG4IHCalDetector::ExtractLayerTowerId(const unsigned int isector, G4VPhysicalVolume *volume)
0266 {
0267 boost::char_separator<char> sep("_");
0268 boost::tokenizer<boost::char_separator<char>> tok(volume->GetName(), sep);
0269 boost::tokenizer<boost::char_separator<char>>::const_iterator tokeniter;
0270 int layer_id = -1, tower_id = -1;
0271 for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
0272 {
0273 if (*tokeniter == "impr")
0274 {
0275 ++tokeniter;
0276 if (tokeniter != tok.end())
0277 {
0278 layer_id = boost::lexical_cast<int>(*tokeniter) / 2;
0279 layer_id--;
0280 if (layer_id < 0 || layer_id >= m_NumScintiPlates)
0281 {
0282 std::cout << "invalid scintillator row " << layer_id
0283 << ", valid range 0 < row < " << m_NumScintiPlates << std::endl;
0284 gSystem->Exit(1);
0285 }
0286 }
0287 else
0288 {
0289 std::cout << PHWHERE << " Error parsing " << volume->GetName()
0290 << " for mother volume number " << std::endl;
0291 gSystem->Exit(1);
0292 }
0293 break;
0294 }
0295 }
0296 for (tokeniter = tok.begin(); tokeniter != tok.end(); ++tokeniter)
0297 {
0298 if (*tokeniter == "pv")
0299 {
0300 ++tokeniter;
0301 if (tokeniter != tok.end())
0302 {
0303 tower_id = boost::lexical_cast<int>(*tokeniter);
0304 }
0305 }
0306 }
0307 int column = map_towerid(tower_id);
0308 int row = map_layerid(layer_id);
0309 return std::make_tuple(isector, row, column);
0310 }
0311
0312
0313 int PHG4IHCalDetector::map_towerid(const int tower_id)
0314 {
0315
0316
0317
0318 int itwr = -1;
0319 int itmp = tower_id / 2;
0320 if (tower_id % 2)
0321 {
0322 itwr = 12 + itmp;
0323 }
0324 else
0325 {
0326 itwr = 11 - itmp;
0327 }
0328 return itwr;
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432 }
0433
0434 int PHG4IHCalDetector::map_layerid(const int layer_id)
0435 {
0436 int rowid = -1;
0437
0438 if (layer_id < 188)
0439 {
0440 rowid = layer_id + 68;
0441 }
0442 else
0443 {
0444 rowid = layer_id - 188;
0445 }
0446
0447 rowid += 4;
0448 if (rowid > 255)
0449 {
0450 rowid -= 256;
0451 }
0452
0453 if (rowid > 255 || rowid < 0)
0454 {
0455 std::cout << PHWHERE << " row id out of range: " << rowid << std::endl;
0456 gSystem->Exit(1);
0457 }
0458 return rowid;
0459 }
0460
0461
0462 void PHG4IHCalDetector::AddGeometryNode()
0463 {
0464 PHNodeIterator iter(topNode());
0465 PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0466 if (!runNode)
0467 {
0468 std::cout << PHWHERE << "Run Node missing, exiting." << std::endl;
0469 gSystem->Exit(1);
0470 exit(1);
0471 }
0472 PHNodeIterator runIter(runNode);
0473 PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_SuperDetector));
0474 if (!RunDetNode)
0475 {
0476 RunDetNode = new PHCompositeNode(m_SuperDetector);
0477 runNode->addNode(RunDetNode);
0478 }
0479 m_TowerGeomNodeName = "TOWERGEOM_" + m_SuperDetector;
0480 m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode(), m_TowerGeomNodeName);
0481 if (!m_RawTowerGeom)
0482 {
0483 m_RawTowerGeom = new RawTowerGeomContainer_Cylinderv1(RawTowerDefs::convert_name_to_caloid(m_SuperDetector));
0484 PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(m_RawTowerGeom, m_TowerGeomNodeName, "PHObject");
0485 RunDetNode->addNode(newNode);
0486 }
0487 double innerrad = m_Params->get_double_param(PHG4HcalDefs::innerrad);
0488 double thickness = m_Params->get_double_param(PHG4HcalDefs::outerrad) - innerrad;
0489 m_RawTowerGeom->set_radius(innerrad);
0490 m_RawTowerGeom->set_thickness(thickness);
0491 m_RawTowerGeom->set_phibins(m_Params->get_int_param(PHG4HcalDefs::n_towers));
0492 m_RawTowerGeom->set_etabins(m_Params->get_int_param("etabins"));
0493 double geom_ref_radius = innerrad + thickness / 2.;
0494 double phistart = m_Params->get_double_param("phistart");
0495 if (!std::isfinite(phistart))
0496 {
0497 std::cout << PHWHERE << " phistart is not finite: " << phistart
0498 << ", exiting now (this will crash anyway)" << std::endl;
0499 gSystem->Exit(1);
0500 }
0501 for (int i = 0; i < m_Params->get_int_param(PHG4HcalDefs::n_towers); i++)
0502 {
0503 double phiend = phistart + 2. * M_PI / m_Params->get_int_param(PHG4HcalDefs::n_towers);
0504 std::pair<double, double> range = std::make_pair(phiend, phistart);
0505 phistart = phiend;
0506 int tempi = i + 1;
0507 if (tempi >= m_Params->get_int_param(PHG4HcalDefs::n_towers))
0508 {
0509 tempi -= m_Params->get_int_param(PHG4HcalDefs::n_towers);
0510 }
0511 m_RawTowerGeom->set_phibounds(tempi, range);
0512 }
0513 double etalowbound = -m_Params->get_double_param("scinti_eta_coverage_neg");
0514 for (int i = 0; i < m_Params->get_int_param("etabins"); i++)
0515 {
0516
0517 double etahibound = etalowbound +
0518 (m_Params->get_double_param("scinti_eta_coverage_neg") + m_Params->get_double_param("scinti_eta_coverage_pos")) / m_Params->get_int_param("etabins");
0519 std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
0520 m_RawTowerGeom->set_etabounds(i, range);
0521 etalowbound = etahibound;
0522 }
0523 for (int iphi = 0; iphi < m_RawTowerGeom->get_phibins(); iphi++)
0524 {
0525 for (int ieta = 0; ieta < m_RawTowerGeom->get_etabins(); ieta++)
0526 {
0527 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::convert_name_to_caloid(m_SuperDetector), ieta, iphi);
0528
0529 const double x(geom_ref_radius * std::cos(m_RawTowerGeom->get_phicenter(iphi)));
0530 const double y(geom_ref_radius * std::sin(m_RawTowerGeom->get_phicenter(iphi)));
0531 const double z(geom_ref_radius / std::tan(PHG4Utils::get_theta(m_RawTowerGeom->get_etacenter(ieta))));
0532
0533 RawTowerGeom *tg = m_RawTowerGeom->get_tower_geometry(key);
0534 if (tg)
0535 {
0536 if (Verbosity() > 0)
0537 {
0538 std::cout << "IHCalDetector::InitRun - Tower geometry " << key << " already exists" << std::endl;
0539 }
0540
0541 if (std::fabs(tg->get_center_x() - x) > 1e-4)
0542 {
0543 std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0544 << std::endl;
0545
0546 return;
0547 }
0548 if (std::fabs(tg->get_center_y() - y) > 1e-4)
0549 {
0550 std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0551 << std::endl;
0552 return;
0553 }
0554 if (std::fabs(tg->get_center_z() - z) > 1e-4)
0555 {
0556 std::cout << "IHCalDetector::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0557 << std::endl;
0558 return;
0559 }
0560 }
0561 else
0562 {
0563 if (Verbosity() > 0)
0564 {
0565 std::cout << "IHCalDetector::InitRun - building tower geometry " << key << "" << std::endl;
0566 }
0567
0568 tg = new RawTowerGeomv1(key);
0569
0570 tg->set_center_x(x);
0571 tg->set_center_y(y);
0572 tg->set_center_z(z);
0573 m_RawTowerGeom->add_tower_geometry(tg);
0574 }
0575 }
0576 }
0577 if (Verbosity() > 0)
0578 {
0579 m_RawTowerGeom->identify();
0580 }
0581 }
