packages/PHField/PHFieldUtility.cc

0001 #include "PHFieldUtility.h"
0002
0003 #include "PHField.h"
0004 #include "PHField2D.h"
0005 #include "PHField3DCartesian.h"
0006 #include "PHField3DCylindrical.h"
0007 #include "PHFieldInterpolated.h"
0008 #include "PHFieldConfig.h"
0009 #include "PHFieldConfigv1.h"
0010 #include "PHFieldUniform.h"
0011
0012 #include <fun4all/Fun4AllServer.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHDataNode.h>
0016 #include <phool/PHIODataNode.h>
0017 #include <phool/PHNodeIterator.h>
0018 #include <phool/PHObject.h>
0019 #include <phool/getClass.h>
0020 #include <phool/phool.h>  // for PHWHERE
0021
0022 #include <TSystem.h>
0023
0024 #include <cassert>
0025 #include <cstdlib>  // for getenv
0026 #include <iostream>
0027
0028 PHField *
0029 PHFieldUtility::BuildFieldMap(const PHFieldConfig *field_config, float inner_radius, float outer_radius, float size_z, const int verbosity)
0030 {
0031   assert(field_config);
0032
0033   if (verbosity)
0034   {
0035     std::cout << "PHFieldUtility::BuildFieldMap - construction field with configuration: ";
0036     field_config->identify();
0037   }
0038
0039   PHField *field(nullptr);
0040
0041   switch (field_config->get_field_config())
0042   {
0043   case PHFieldConfig::kFieldUniform:
0044     //    return "Constant field";
0045
0046     field = new PHFieldUniform(
0047         field_config->get_field_mag_x(),
0048         field_config->get_field_mag_y(),
0049         field_config->get_field_mag_z());
0050
0051     break;
0052   case PHFieldConfig::kField2D:
0053     //    return "2D field map expressed in cylindrical coordinates";
0054     field = new PHField2D(
0055         field_config->get_filename(),
0056         verbosity,
0057         field_config->get_magfield_rescale());
0058     break;
0059
0060   case PHFieldConfig::kField3DCylindrical:
0061     //    return "3D field map expressed in cylindrical coordinates";
0062     field = new PHField3DCylindrical(
0063         field_config->get_filename(),
0064         verbosity,
0065         field_config->get_magfield_rescale());
0066     break;
0067
0068   case PHFieldConfig::Field3DCartesian:
0069     //    return "3D field map expressed in Cartesian coordinates";
0070     field = new PHField3DCartesian(
0071         field_config->get_filename(),
0072         field_config->get_magfield_rescale(),
0073         inner_radius,
0074         outer_radius,
0075         size_z);
0076     break;
0077   case PHFieldConfig::FieldInterpolated:
0078     //    return "3d interpolated fieldmap"
0079     field = new PHFieldInterpolated;
0080     field->Verbosity(1); // Will print to std::cout when loading
0081     try {
0082       dynamic_cast<PHFieldInterpolated*>(field)->load_fieldmap (
0083         field_config->get_filename(),
0084         field_config->get_magfield_rescale()
0085       );
0086     } catch (std::exception const& e) {
0087       delete field;
0088       field = nullptr;
0089       std::cout << e.what() << std::endl;
0090       gSystem->Exit(1);
0091     }
0092     break;
0093
0094   default:
0095     std::cout << "PHFieldUtility::BuildFieldMap - Invalid Field Configuration: " << field_config->get_field_config() << std::endl;
0096     gSystem->Exit(1);
0097   }
0098   assert(field);  // Check for Invalid Field
0099   return field;
0100 }
0101
0102 //! Make a default PHFieldConfig
0103 //! Field map = /phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root
0104 //! Field Scale to 1.4/1.5
0105 //! \output default field configuration object. Caller assumes ownership
0106 PHFieldConfig *PHFieldUtility::DefaultFieldConfig()
0107 {
0108   char *calibrationroot = getenv("CALIBRATIONROOT");
0109   std::string fieldmap = "sphenix3dbigmapxyz_gap_rebuild.root";
0110   if (calibrationroot != nullptr)
0111   {
0112     fieldmap = std::string(calibrationroot) + "/Field/Map/" + fieldmap;
0113   }
0114   return new PHFieldConfigv1(PHFieldConfigv1::Field3DCartesian, fieldmap, 1.);
0115 }
0116
0117 //! Get transient PHField from DST nodes. If not found, make a new one based on default_config
0118 PHField *
0119 PHFieldUtility::GetFieldMapNode(const PHFieldConfig *default_config, PHCompositeNode *topNode, const int verbosity)
0120 {
0121   if (topNode == nullptr)
0122   {
0123     topNode = Fun4AllServer::instance()->topNode();
0124   }
0125   assert(topNode);
0126   PHNodeIterator iter(topNode);
0127
0128   // Looking for the RUN node
0129   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0130   if (!parNode)
0131   {
0132     std::cout << PHWHERE << ": PAR Node missing, request aborting.";
0133     gSystem->Exit(1);
0134   }
0135
0136   PHField *field = findNode::getClass<PHField>(parNode, GetDSTFieldMapNodeName());
0137   if (!field)
0138   {
0139     PHFieldConfig *field_config = GetFieldConfigNode(default_config, topNode, verbosity);
0140     assert(field_config);
0141
0142     field = BuildFieldMap(field_config, verbosity > 0 ? verbosity - 1 : verbosity);
0143     assert(field);
0144
0145     parNode->addNode(new PHDataNode<PHField>(field, GetDSTFieldMapNodeName()));
0146   }
0147
0148   return field;
0149 }
0150
0151 //! Get persistent PHGeomIOTGeo from DST nodes. If not found, make a new one
0152 PHFieldConfig *
0153 PHFieldUtility::GetFieldConfigNode(const PHFieldConfig *default_config, PHCompositeNode *topNode, const int verbosity)
0154 {
0155   if (topNode == nullptr)
0156   {
0157     topNode = Fun4AllServer::instance()->topNode();
0158   }
0159   assert(topNode);
0160
0161   PHNodeIterator iter(topNode);
0162
0163   // Looking for the RUN node
0164   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0165   if (!runNode)
0166   {
0167     std::cout << PHWHERE << ": RUN Node missing, aborting.";
0168     gSystem->Exit(1);
0169   }
0170
0171   PHFieldConfig *field = findNode::getClass<PHFieldConfig>(runNode, GetDSTConfigNodeName());
0172   if (!field)
0173   {
0174     if (!default_config)
0175     {
0176       field = DefaultFieldConfig();
0177       if (verbosity)
0178       {
0179         std::cout << "PHFieldUtility::GetFieldConfigNode - field map with configuration from build-in default: ";
0180         field->identify();
0181       }
0182     }
0183     else
0184     {
0185       field = dynamic_cast<PHFieldConfig *>(default_config->CloneMe());
0186       if (verbosity)
0187       {
0188         std::cout << "PHFieldUtility::GetFieldConfigNode - field map with configuration from input default: ";
0189         field->identify();
0190       }
0191     }
0192
0193     assert(field);
0194     runNode->addNode(new PHIODataNode<PHObject>(field, GetDSTConfigNodeName(), "PHObject"));
0195   }
0196   else
0197   {
0198     if (verbosity)
0199     {
0200       std::cout << "PHFieldUtility::GetFieldConfigNode - field map with configuration from DST/RUN: ";
0201       field->identify();
0202     }
0203   }
0204
0205   return field;
0206 }