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 using namespace std;
0030
0031 //! DST node -> TGeoManager for downstream use
0032 TGeoManager *
0033 PHGeomUtility::GetTGeoManager(PHCompositeNode *topNode)
0034 {
0035   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode, true);
0036   if (!dst_geom)
0037   {
0038     cout << __PRETTY_FUNCTION__
0039          << " - Error - Can NOT construct geometry node." << endl;
0040     exit(1);
0041     return nullptr;
0042   }
0043
0044   if (not dst_geom->isValid())
0045   {
0046     // try to construct the geometry node
0047     dst_geom = LoadFromIONode(topNode);
0048   }
0049
0050   if (TGeoManager::GetDefaultUnits() != TGeoManager::kRootUnits )
0051   {
0052     cout << __PRETTY_FUNCTION__ << " TGeoManager was not constructed with RootUnits, which potentially leads to unit mismatch with Fun4All. This is considered a fatal error."
0053         <<endl;
0054
0055     exit(1);
0056     return nullptr;
0057   }
0058
0059   UpdateIONode(topNode);
0060
0061   return dst_geom->GetGeometry();
0062 }
0063
0064 void PHGeomUtility::ExportGeomtry(PHCompositeNode *topNode,
0065                                   const std::string &geometry_file)
0066 {
0067   TGeoManager *tgeo = GetTGeoManager(topNode);
0068
0069   assert(tgeo);
0070
0071   tgeo->Export(geometry_file.c_str());
0072 }
0073
0074 int PHGeomUtility::ImportGeomFile(PHCompositeNode *topNode,
0075                                   const std::string &geometry_file)
0076 {
0077   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode);
0078   assert(dst_geom);
0079
0080   dst_geom->Reset();
0081
0082   TGeoManager::SetVerboseLevel(GetVerbosity());
0083
0084   // force TGeoManager to use the Fun4All unit of cm
0085 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,23,2)
0086   TGeoManager::LockDefaultUnits(kFALSE);
0087   TGeoManager::SetDefaultUnits( TGeoManager::kRootUnits );
0088   TGeoManager::LockDefaultUnits(kTRUE);
0089 #else
0090   TGeoManager::SetDefaultRootUnits();
0091 #endif
0092
0093   dst_geom->SetGeometry(TGeoManager::Import(geometry_file.c_str()));
0094
0095   if (dst_geom->GetGeometry() == nullptr)
0096   {
0097     cout << __PRETTY_FUNCTION__ << "failed to import " << geometry_file
0098          << endl;
0099     return Fun4AllReturnCodes::ABORTRUN;
0100   }
0101
0102   UpdateIONode(topNode);
0103
0104   return Fun4AllReturnCodes::EVENT_OK;
0105 }
0106
0107 int PHGeomUtility::ImportCurrentTGeoManager(PHCompositeNode *topNode)
0108 {
0109   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode);
0110   assert(dst_geom);
0111
0112   if (dst_geom->GetGeometry() == gGeoManager)
0113     return Fun4AllReturnCodes::EVENT_OK;  // noting to be done
0114
0115   assert(not 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     ostringstream serr;
0134     serr << __PRETTY_FUNCTION__ << ": PAR Node missing, request aborting.";
0135     cout << serr.str() << endl;
0136
0137     throw runtime_error(serr.str());
0138
0139     return nullptr;
0140   }
0141
0142   PHGeomTGeo *dst_geom = findNode::getClass<PHGeomTGeo>(parNode,
0143                                                         GetDSTNodeName());
0144   if (!dst_geom and 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     ostringstream serr;
0167     serr << __PRETTY_FUNCTION__ << ": RUN Node missing, request aborting.";
0168     cout << serr.str() << endl;
0169
0170     throw runtime_error(serr.str());
0171
0172     return nullptr;
0173   }
0174
0175   PHGeomIOTGeo *dst_geom = findNode::getClass<PHGeomIOTGeo>(runNode,
0176                                                             GetDSTIONodeName());
0177   if (!dst_geom and 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   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   fstream ifile(file_name, ios_base::in);
0211
0212   if (ifile)
0213   {
0214     ifile.close();
0215     if (remove(file_name.c_str()) != 0)
0216     {
0217       cout << __PRETTY_FUNCTION__ << " - Error - can not remove file "
0218            << file_name << endl;
0219       return false;
0220     }
0221     else
0222       return true;
0223   }
0224   else
0225     return true;  // file do not exist
0226 }
0227
0228 //! Update persistent PHGeomIOTGeo node RUN/GEOMETRY_IO based on run-time object PHGeomTGeo at RUN/GEOMETRY
0229 //! \return the updated PHGeomIOTGeo from DST tree
0230 PHGeomIOTGeo *
0231 PHGeomUtility::UpdateIONode(PHCompositeNode *topNode)
0232 {
0233   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode, false);
0234
0235   if (not dst_geom)
0236   {
0237     cout << __PRETTY_FUNCTION__
0238          << " - ERROR - failed to update PHGeomIOTGeo node RUN/GEOMETRY_IO due to missing PHGeomTGeo node at RUN/GEOMETRY"
0239          << endl;
0240     return nullptr;
0241   }
0242   if (not dst_geom->isValid())
0243   {
0244     cout << __PRETTY_FUNCTION__
0245          << " - ERROR - failed to update PHGeomIOTGeo node RUN/GEOMETRY_IO due to invalid PHGeomTGeo node at RUN/GEOMETRY"
0246          << endl;
0247     return nullptr;
0248   }
0249
0250   PHGeomIOTGeo *dst_geom_io = GetGeomIOTGeoNode(topNode, true);
0251   assert(dst_geom_io);
0252
0253   dst_geom_io->SetGeometry(dst_geom->GetGeometry()->GetTopVolume());
0254
0255   return dst_geom_io;
0256 }
0257
0258 //! Build or update PHGeomTGeo node RUN/GEOMETRY based on the persistent PHGeomIOTGeo node RUN/GEOMETRY_IO
0259 //! \return the updated PHGeomTGeo from DST tree
0260 PHGeomTGeo *
0261 PHGeomUtility::LoadFromIONode(PHCompositeNode *topNode)
0262 {
0263   PHGeomIOTGeo *dst_geom_io = GetGeomIOTGeoNode(topNode, false);
0264
0265   if (not dst_geom_io)
0266   {
0267     cout << __PRETTY_FUNCTION__
0268          << " - ERROR - failed to update PHGeomTGeo node RUN/GEOMETRY due to missing PHGeomIOTGeo node at RUN/GEOMETRY_IO"
0269          << endl;
0270     return nullptr;
0271   }
0272   if (not dst_geom_io->isValid())
0273   {
0274     cout << __PRETTY_FUNCTION__
0275          << " - ERROR - failed to update PHGeomTGeo node RUN/GEOMETRY due to invalid PHGeomIOTGeo node at RUN/GEOMETRY_IO"
0276          << endl;
0277     return nullptr;
0278   }
0279
0280   // build new TGeoManager
0281   TGeoManager::SetVerboseLevel(GetVerbosity());
0282   TGeoManager *tgeo = dst_geom_io->ConstructTGeoManager();
0283   tgeo->CloseGeometry();
0284
0285   PHGeomTGeo *dst_geom = GetGeomTGeoNode(topNode, true);
0286   assert(dst_geom);
0287   dst_geom->SetGeometry(tgeo);
0288
0289   return dst_geom;
0290 }
0291
0292 //! Verbosity for geometry IO like, TGeoMangers
0293 void PHGeomUtility::SetVerbosity(int v)
0294 {
0295   recoConsts *rc = recoConsts::instance();
0296   rc->set_IntFlag("PHGEOMETRY_VERBOSITY", v);
0297 }
0298
0299 //! Verbosity for geometry IO like, TGeoMangers
0300 int PHGeomUtility::GetVerbosity()
0301 {
0302   recoConsts *rc = recoConsts::instance();
0303   if (rc->FlagExist("PHGEOMETRY_VERBOSITY"))
0304     return rc->get_IntFlag("PHGEOMETRY_VERBOSITY");
0305   else
0306     return 0;
0307 }