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