packages/PHGeometry/PHGeomUtility.cc

0001 #include "PHGeomUtility.h"
0002
0003 #include "PHGeomIOTGeo.h"
0004 #include "PHGeomTGeo.h"
0005
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/PHDataNode.h>
0010 #include <phool/PHIODataNode.h>
0011 #include <phool/PHNodeIterator.h>
0012 #include <phool/PHObject.h>
0013 #include <phool/getClass.h>
0014 #include <phool/recoConsts.h>
0015
0016 #include <TGeoManager.h>
0017
0018 #include <uuid/uuid.h>
0019
0020 #include <unistd.h>  // for generate unique local file
0021 #include <cassert>
0022 #include <cstdio>
0023 #include <cstdlib>
0024 #include <fstream>
0025 #include <iostream>
0026 #include <sstream>
0027 #include <stdexcept>
0028
0029 //! DST node -> TGeoManager for downstream use
0030 TGeoManager *
0031 PHGeomUtility::GetTGeoManager(PHCompositeNode *topNode)
0032 {
0033   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode, true);
0034   if (!dst_geom)
0035   {
0036     std::cout << __PRETTY_FUNCTION__
0037               << " - Error - Can NOT construct geometry node." << std::endl;
0038     exit(1);
0039     return nullptr;
0040   }
0041
0042   if (!dst_geom->isValid())
0043   {
0044     // try to construct the geometry node
0045     dst_geom = LoadFromIONode(topNode);
0046   }
0047
0048   if (TGeoManager::GetDefaultUnits() != TGeoManager::kRootUnits)
0049   {
0050     std::cout << __PRETTY_FUNCTION__ << " TGeoManager was not constructed with RootUnits, which potentially leads to unit mismatch with Fun4All. This is considered a fatal error."
0051               << std::endl;
0052
0053     exit(1);
0054     return nullptr;
0055   }
0056
0057   UpdateIONode(topNode);
0058
0059   return dst_geom->GetGeometry();
0060 }
0061
0062 void PHGeomUtility::ExportGeomtry(PHCompositeNode *topNode,
0063                                   const std::string &geometry_file)
0064 {
0065   TGeoManager *tgeo = GetTGeoManager(topNode);
0066
0067   assert(tgeo);
0068
0069   tgeo->Export(geometry_file.c_str());
0070 }
0071
0072 int PHGeomUtility::ImportGeomFile(PHCompositeNode *topNode,
0073                                   const std::string &geometry_file)
0074 {
0075   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode);
0076   assert(dst_geom);
0077
0078   dst_geom->Reset();
0079
0080   TGeoManager::SetVerboseLevel(GetVerbosity());
0081
0082   // force TGeoManager to use the Fun4All unit of cm
0083 #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 23, 2)
0084   TGeoManager::LockDefaultUnits(kFALSE);
0085   TGeoManager::SetDefaultUnits(TGeoManager::kRootUnits);
0086   TGeoManager::LockDefaultUnits(kTRUE);
0087 #else
0088   TGeoManager::SetDefaultRootUnits();
0089 #endif
0090
0091   dst_geom->SetGeometry(TGeoManager::Import(geometry_file.c_str()));
0092
0093   if (dst_geom->GetGeometry() == nullptr)
0094   {
0095     std::cout << __PRETTY_FUNCTION__ << "failed to import " << geometry_file
0096               << std::endl;
0097     return Fun4AllReturnCodes::ABORTRUN;
0098   }
0099
0100   UpdateIONode(topNode);
0101
0102   return Fun4AllReturnCodes::EVENT_OK;
0103 }
0104
0105 int PHGeomUtility::ImportCurrentTGeoManager(PHCompositeNode *topNode)
0106 {
0107   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode);
0108   assert(dst_geom);
0109
0110   if (dst_geom->GetGeometry() == gGeoManager)
0111   {
0112     return Fun4AllReturnCodes::EVENT_OK;  // noting to be done
0113   }
0114
0115   assert(!dst_geom->isValid());  // check that it is uninitialized
0116   dst_geom->SetGeometry(gGeoManager);
0117   TGeoManager::SetVerboseLevel(GetVerbosity());
0118
0119   return Fun4AllReturnCodes::EVENT_OK;
0120 }
0121
0122 //! Get non-persistent PHGeomTGeo from DST nodes. If not found, make a new one
0123 PHGeomTGeo *
0124 PHGeomUtility::GetGeomTGeoNode(PHCompositeNode *topNode, bool build_new)
0125 {
0126   PHNodeIterator iter(topNode);
0127
0128   // Looking for the RUN node
0129   PHCompositeNode *parNode = static_cast<PHCompositeNode *>(iter.findFirst(
0130       "PHCompositeNode", "PAR"));
0131   if (!parNode)
0132   {
0133     std::ostringstream serr;
0134     serr << __PRETTY_FUNCTION__ << ": PAR Node missing, request aborting.";
0135     std::cout << serr.str() << std::endl;
0136
0137     throw std::runtime_error(serr.str());
0138
0139     return nullptr;
0140   }
0141
0142   PHGeomTGeo *dst_geom = findNode::getClass<PHGeomTGeo>(parNode,
0143                                                         GetDSTNodeName());
0144   if (!dst_geom && build_new)
0145   {
0146     dst_geom = new PHGeomTGeo();
0147     PHDataNode<PHObject> *GeomNode = new PHDataNode<PHObject>(dst_geom,
0148                                                               GetDSTNodeName(), "PHObject");
0149     parNode->addNode(GeomNode);
0150   }
0151
0152   return dst_geom;
0153 }
0154
0155 //! Get persistent PHGeomIOTGeo from DST nodes. If not found, make a new one
0156 PHGeomIOTGeo *
0157 PHGeomUtility::GetGeomIOTGeoNode(PHCompositeNode *topNode, bool build_new)
0158 {
0159   PHNodeIterator iter(topNode);
0160
0161   // Looking for the RUN node
0162   PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst(
0163       "PHCompositeNode", "RUN"));
0164   if (!runNode)
0165   {
0166     std::ostringstream serr;
0167     serr << __PRETTY_FUNCTION__ << ": RUN Node missing, request aborting.";
0168     std::cout << serr.str() << std::endl;
0169
0170     throw std::runtime_error(serr.str());
0171
0172     return nullptr;
0173   }
0174
0175   PHGeomIOTGeo *dst_geom = findNode::getClass<PHGeomIOTGeo>(runNode,
0176                                                             GetDSTIONodeName());
0177   if (!dst_geom && build_new)
0178   {
0179     dst_geom = new PHGeomIOTGeo();
0180     PHIODataNode<PHObject> *GeomNode = new PHIODataNode<PHObject>(dst_geom,
0181                                                                   GetDSTIONodeName(), "PHObject");
0182     runNode->addNode(GeomNode);
0183   }
0184
0185   return dst_geom;
0186 }
0187
0188 std::string
0189 PHGeomUtility::GenerateGeometryFileName(const std::string &filename_extension)
0190 {
0191   std::ostringstream file;
0192
0193   uuid_t uu;
0194   uuid_generate(uu);
0195   char uuid[50];
0196   uuid_unparse(uu, uuid);
0197
0198   file << mg_GenerateGeometryFileNameBase << "/"
0199        << "PHGeomUtility_geom_file_" << uuid << "."
0200        << filename_extension;
0201
0202   return file.str();
0203 }
0204
0205 std::string PHGeomUtility::mg_GenerateGeometryFileNameBase = "/tmp";
0206
0207 //! delete the geometry file after use
0208 bool PHGeomUtility::RemoveGeometryFile(const std::string &file_name)
0209 {
0210   std::fstream ifile(file_name, std::ios_base::in);
0211
0212   if (ifile)
0213   {
0214     ifile.close();
0215     if (remove(file_name.c_str()) != 0)
0216     {
0217       std::cout << __PRETTY_FUNCTION__ << " - Error - can not remove file "
0218                 << file_name << std::endl;
0219       return false;
0220     }
0221     return true;
0222   }
0223   return true;  // file do not exist
0224 }
0225
0226 //! Update persistent PHGeomIOTGeo node RUN/GEOMETRY_IO based on run-time object PHGeomTGeo at RUN/GEOMETRY
0227 //! \return the updated PHGeomIOTGeo from DST tree
0228 PHGeomIOTGeo *
0229 PHGeomUtility::UpdateIONode(PHCompositeNode *topNode)
0230 {
0231   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode, false);
0232
0233   if (!dst_geom)
0234   {
0235     std::cout << __PRETTY_FUNCTION__
0236               << " - ERROR - failed to update PHGeomIOTGeo node RUN/GEOMETRY_IO due to missing PHGeomTGeo node at RUN/GEOMETRY"
0237               << std::endl;
0238     return nullptr;
0239   }
0240   if (!dst_geom->isValid())
0241   {
0242     std::cout << __PRETTY_FUNCTION__
0243               << " - ERROR - failed to update PHGeomIOTGeo node RUN/GEOMETRY_IO due to invalid PHGeomTGeo node at RUN/GEOMETRY"
0244               << std::endl;
0245     return nullptr;
0246   }
0247
0248   PHGeomIOTGeo *dst_geom_io = GetGeomIOTGeoNode(topNode, true);
0249   assert(dst_geom_io);
0250
0251   dst_geom_io->SetGeometry(dst_geom->GetGeometry()->GetTopVolume());
0252
0253   return dst_geom_io;
0254 }
0255
0256 //! Build or update PHGeomTGeo node RUN/GEOMETRY based on the persistent PHGeomIOTGeo node RUN/GEOMETRY_IO
0257 //! \return the updated PHGeomTGeo from DST tree
0258 PHGeomTGeo *
0259 PHGeomUtility::LoadFromIONode(PHCompositeNode *topNode)
0260 {
0261   PHGeomIOTGeo *dst_geom_io = GetGeomIOTGeoNode(topNode, false);
0262
0263   if (!dst_geom_io)
0264   {
0265     std::cout << __PRETTY_FUNCTION__
0266               << " - ERROR - failed to update PHGeomTGeo node RUN/GEOMETRY due to missing PHGeomIOTGeo node at RUN/GEOMETRY_IO"
0267               << std::endl;
0268     return nullptr;
0269   }
0270   if (!dst_geom_io->isValid())
0271   {
0272     std::cout << __PRETTY_FUNCTION__
0273               << " - ERROR - failed to update PHGeomTGeo node RUN/GEOMETRY due to invalid PHGeomIOTGeo node at RUN/GEOMETRY_IO"
0274               << std::endl;
0275     return nullptr;
0276   }
0277
0278   // build new TGeoManager
0279   TGeoManager::SetVerboseLevel(GetVerbosity());
0280   TGeoManager *tgeo = dst_geom_io->ConstructTGeoManager();
0281   tgeo->CloseGeometry();
0282
0283   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode, true);
0284   assert(dst_geom);
0285   dst_geom->SetGeometry(tgeo);
0286
0287   return dst_geom;
0288 }
0289
0290 //! Verbosity for geometry IO like, TGeoMangers
0291 void PHGeomUtility::SetVerbosity(int v)
0292 {
0293   recoConsts *rc = recoConsts::instance();
0294   rc->set_IntFlag("PHGEOMETRY_VERBOSITY", v);
0295 }
0296
0297 //! Verbosity for geometry IO like, TGeoMangers
0298 int PHGeomUtility::GetVerbosity()
0299 {
0300   recoConsts *rc = recoConsts::instance();
0301   if (rc->FlagExist("PHGEOMETRY_VERBOSITY"))
0302   {
0303     return rc->get_IntFlag("PHGEOMETRY_VERBOSITY");
0304   }
0305   return 0;
0306 }