packages/trackreco/MakeActsGeometry.h

0001 /*!
0002  *  \file       MakeActsGeometry.h
0003  *  \brief      Make Acts geometry from sPHENIX TGeometry
0004  *  \details    Make Acts geometry from sPHENIX TGeometry
0005  *  \author     Tony Frawley <afrawley@fsu.edu>
0006  */
0007
0008 #ifndef MAKE_ACTS_GEOMETRY_H
0009 #define MAKE_ACTS_GEOMETRY_H
0010
0011 #include <fun4all/SubsysReco.h>
0012 #include <trackbase/TrkrDefs.h>
0013
0014 #include <trackbase/ActsGeometry.h>
0015
0016 #include <Acts/Definitions/Algebra.hpp>
0017 #include <Acts/EventData/MeasurementHelpers.hpp>
0018 #include <Acts/Geometry/TrackingGeometry.hpp>
0019 #include <Acts/MagneticField/MagneticFieldContext.hpp>
0020 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0021 #include <Acts/Utilities/BinnedArray.hpp>
0022 #include <Acts/Utilities/CalibrationContext.hpp>
0023 #include <Acts/Utilities/Logger.hpp>
0024
0025 #include <trackbase/TGeoDetectorWithOptions.h>
0026 #ifndef __CLING__
0027 #include <boost/program_options.hpp>
0028 #endif
0029 #include <map>
0030 #include <memory>
0031 #include <string>
0032 #include <vector>
0033
0034 class PHCompositeNode;
0035 class PHG4CylinderGeomContainer;
0036 class PHG4TpcCylinderGeomContainer;
0037 class TGeoManager;
0038 class TGeoNode;
0039 class TGeoVolume;
0040
0041 namespace Acts
0042 {
0043   class Surface;
0044 }
0045
0046 using Surface = std::shared_ptr<const Acts::Surface>;
0047 using TrackingGeometry = std::shared_ptr<const Acts::TrackingGeometry>;
0048 // using TrackingGeometry = std::shared_ptr<Acts::TrackingGeometry>;
0049 using TrackingVolumePtr = std::shared_ptr<const Acts::TrackingVolume>;
0050
0051 /**
0052  * This class is responsible for building the ActsGeometry from the sPHENIX
0053  * TGeometry. The code largely follows examples within the ACTFW code,
0054  * specifically GeometryExampleBase.cpp.
0055  * This class puts several nodes on the node tree which relate Acts::Surfaces
0056  * to sPHENIX TGeo objects, for use in building Acts SourceLinks.
0057  */
0058 class MakeActsGeometry : public SubsysReco
0059 {
0060  public:
0061   //! Default constructor
0062   MakeActsGeometry(const std::string &name = "MakeActsGeometry");
0063
0064   //! Destructor
0065   ~MakeActsGeometry() override = default;
0066
0067   int Init(PHCompositeNode *topNode) override;
0068   int InitRun(PHCompositeNode *topNode) override;
0069
0070   void loadMagField(const bool field) { m_useField = field; }
0071   void setMagField(const std::string &magField)
0072   {
0073     m_magField = magField;
0074   }
0075   void setMagFieldRescale(double magFieldRescale)
0076   {
0077     m_magFieldRescale = magFieldRescale;
0078   }
0079
0080   // void useInttSurveyGeom(const bool useSurveyGeom) { m_useInttSurveyGeom = useSurveyGeom; }
0081
0082   void setMvtxDev(double array[6])
0083   {
0084     m_mvtxDevs[0] = array[0];
0085     m_mvtxDevs[1] = array[1];
0086     m_mvtxDevs[2] = array[2];
0087     m_mvtxDevs[3] = array[3];
0088     m_mvtxDevs[4] = array[4];
0089     m_mvtxDevs[5] = array[5];
0090
0091     mvtxParam = true;
0092   }
0093   void setInttDev(double array[6])
0094   {
0095     m_inttDevs[0] = array[0];
0096     m_inttDevs[1] = array[1];
0097     m_inttDevs[2] = array[2];
0098     m_inttDevs[3] = array[3];
0099     m_inttDevs[4] = array[4];
0100     m_inttDevs[5] = array[5];
0101
0102     inttParam = true;
0103   }
0104   void setTpcDev(double array[6])
0105   {
0106     m_tpcDevs[0] = array[0];
0107     m_tpcDevs[1] = array[1];
0108     m_tpcDevs[2] = array[2];
0109     m_tpcDevs[3] = array[3];
0110     m_tpcDevs[4] = array[4];
0111     m_tpcDevs[5] = array[5];
0112
0113     tpcParam = true;
0114   }
0115   void setMmDev(double array[6])
0116   {
0117     m_mmDevs[0] = array[0];
0118     m_mmDevs[1] = array[1];
0119     m_mmDevs[2] = array[2];
0120     m_mmDevs[3] = array[3];
0121     m_mmDevs[4] = array[4];
0122     m_mmDevs[5] = array[5];
0123
0124     mmParam = true;
0125   }
0126
0127   void misalignmentFactor(uint8_t layer, const double misalignment)
0128   {
0129     auto it = m_misalignmentFactor.find(layer);
0130     if (it != m_misalignmentFactor.end())
0131     {
0132       it->second = misalignment;
0133       return;
0134     }
0135   }
0136
0137   double getSurfStepPhi() { return m_surfStepPhi; }
0138   double getSurfStepZ() { return m_surfStepZ; }
0139
0140   void set_drift_velocity(double vd) { m_drift_velocity = vd; }
0141   void set_tpc_tzero(double tz) { m_tpc_tzero = tz; }
0142
0143   void set_nSurfPhi(unsigned int value)
0144   {
0145     m_nSurfPhi = value;
0146   }
0147
0148   void set_mvtx_applymisalign(bool b) { m_mvtxapplymisalign = b; }
0149   void set_intt_survey(bool surv) { m_inttSurvey = surv; }
0150
0151  private:
0152   /// Main function to build all acts geometry for use in the fitting modules
0153   int buildAllGeometry(PHCompositeNode *topNode);
0154
0155   //! Get all the nodes
0156   int getNodes(PHCompositeNode *);
0157
0158   //! Create New nodes
0159   int createNodes(PHCompositeNode *);
0160
0161   /// Functions to edit TGeoManager to include TPC boxes
0162   void setPlanarSurfaceDivisions();
0163   void editTPCGeometry(PHCompositeNode *topNode);
0164   void addActsTpcSurfaces(TGeoVolume *tpc_gas_vol,
0165                           TGeoManager *geoManager);
0166
0167   /// Silicon layers made by BuildSiliconLayers and its helper functions
0168   void buildActsSurfaces();
0169
0170   /// Function that mimics ActsExamples::GeometryExampleBase
0171   void makeGeometry(int argc, char *argv[],
0172                     ActsExamples::TGeoDetectorWithOptions &detector);
0173 #ifndef __CLING__
0174   std::pair<std::shared_ptr<const Acts::TrackingGeometry>,
0175             std::vector<std::shared_ptr<ActsExamples::IContextDecorator>>>
0176   build(const boost::program_options::variables_map &vm,
0177         ActsExamples::TGeoDetectorWithOptions &detector);
0178 #endif
0179   void readTGeoLayerBuilderConfigsFile(const std::string &path,
0180                                        ActsExamples::TGeoDetector::Config &config);
0181
0182   void setMaterialResponseFile(std::string &responseFile,
0183                                std::string &materialFile);
0184
0185   /// Get hitsetkey from TGeoNode for each detector geometry
0186   void getInttKeyFromNode(TGeoNode *gnode);
0187   void getMvtxKeyFromNode(TGeoNode *gnode);
0188   void getTpcKeyFromNode(TGeoNode *gnode);
0189
0190   /// Make the Surface<-->TrkrDef::hitsetkey map pairs for each of
0191   /// the various subdetectors
0192   void makeMvtxMapPairs(TrackingVolumePtr &mvtxVolume);
0193   void makeInttMapPairs(TrackingVolumePtr &inttVolume);
0194   void makeTpcMapPairs(TrackingVolumePtr &tpcVolume);
0195
0196   /// bind micromegas surfaces to hitset id
0197   void makeMmMapPairs(TrackingVolumePtr &tpcVolume);
0198
0199   /// Get subdetector hitsetkey from the local sensor unit coordinates
0200   TrkrDefs::hitsetkey getMvtxHitSetKeyFromCoords(unsigned int layer,
0201                                                  std::vector<double> &world);
0202   TrkrDefs::hitsetkey getInttHitSetKeyFromCoords(unsigned int layer,
0203                                                  std::vector<double> &world);
0204   TrkrDefs::hitsetkey getTpcHitSetKeyFromCoords(std::vector<double> &world);
0205
0206   //   /// Helper diagnostic function for identifying active layers in subdetectors
0207   //   void isActive(TGeoNode *gnode, int nmax_print);
0208
0209   //   /// Makes map of TrkrHitSetKey<-->TGeoNode
0210   //   void makeTGeoNodeMap(PHCompositeNode *topNode);
0211
0212   void unpackVolumes();
0213
0214   /// Subdetector geometry containers for getting layer information
0215   PHG4CylinderGeomContainer *m_geomContainerMvtx = nullptr;
0216   PHG4CylinderGeomContainer *m_geomContainerIntt = nullptr;
0217   PHG4CylinderGeomContainer *m_geomContainerMicromegas = nullptr;
0218   PHG4TpcCylinderGeomContainer *m_geomContainerTpc = nullptr;
0219   TGeoManager *m_geoManager = nullptr;
0220
0221   // Switch to use or not use the INTT survey geometry
0222   bool m_inttSurvey = true;
0223   const float m_inttbarrelcenter_survey_x = 0.4026857142857132 / 10.;
0224   const float m_inttbarrelcenter_survey_y = -2.886627321428573 / 10.;
0225
0226   // Switch for applying misalignment for mvtx
0227   bool m_mvtxapplymisalign = false;
0228   std::vector<double> v_globaldisplacement = {0., 0., 0.};
0229
0230   bool m_useField = true;
0231   std::map<uint8_t, double> m_misalignmentFactor;
0232
0233   /// Several maps that connect Acts world to sPHENIX G4 world
0234   std::map<TrkrDefs::hitsetkey, TGeoNode *> m_clusterNodeMap;
0235   std::map<TrkrDefs::hitsetkey, Surface> m_clusterSurfaceMapSilicon;
0236   std::map<unsigned int, std::vector<Surface>> m_clusterSurfaceMapTpcEdit;  // uses layer as key
0237   std::map<TrkrDefs::hitsetkey, Surface> m_clusterSurfaceMapMmEdit;
0238
0239   /// These don't change, we are building the tpc this way!
0240   static constexpr unsigned int m_nTpcLayers = 48;
0241   static constexpr unsigned int m_nTpcModulesPerLayer = 12;
0242   static constexpr unsigned int m_nTpcSides = 2;
0243
0244   /// TPC Acts::Surface subdivisions
0245   double m_minSurfZ = 0.;
0246   /// This value must be less than the TPC gas volume in TGeo, which
0247   /// is 105.22 cm
0248   double m_maxSurfZ = 105.42;
0249   unsigned int m_nSurfZ = 1;
0250   unsigned int m_nSurfPhi = 12;
0251   double m_surfStepPhi = 0;
0252   double m_surfStepZ = 0;
0253   double m_moduleStepPhi = 0;
0254   double m_modulePhiStart = 0;
0255
0256   /// Debugger for printing out tpc active volumes
0257   //  int nprint_tpc = 0;
0258
0259   /// TPC TGeoManager editing box surfaces subdivisions
0260   const static int m_nTpcSectors = 3;
0261
0262   double m_layerRadius[m_nTpcLayers] = {0};
0263   double m_layerThickness[m_nTpcLayers] = {0};
0264
0265   // Spaces to prevent boxes from touching when placed
0266   const double half_width_clearance_thick = 0.4999;
0267   const double half_width_clearance_phi = 0.4999;
0268   /// z does not need spacing as the boxes are rotated around the z axis
0269   const double half_width_clearance_z = 0.5;
0270
0271   /// The acts geometry object
0272   ActsExamples::TGeoDetectorWithOptions m_detector;
0273
0274   /// Acts geometry objects that are needed to create (for example) the fitter
0275   TrackingGeometry m_tGeometry;
0276   std::shared_ptr<Acts::MagneticFieldProvider> m_magneticField;
0277   Acts::GeometryContext m_geoCtxt;
0278
0279   /// Structs to put on the node tree which carry around ActsGeom info
0280   ActsGeometry *m_actsGeometry = nullptr;
0281
0282   std::map<unsigned int, unsigned int> base_layer_map = {{10, 0}, {12, 3}, {14, 7}, {16, 55}};
0283   unsigned int mvtx_chips_per_stave = 9;
0284
0285   /// Verbosity value handed from PHActsSourceLinks
0286   //  int m_verbosity = 0;
0287
0288   double m_drift_velocity = 8.0e-03;  // cm/ns, override from macro
0289   double m_tpc_tzero = 0.0;  // ns, override from macro
0290
0291   /// Magnetic field components to set Acts magnetic field
0292   std::string m_magField = "1.4";
0293   double m_magFieldRescale = -1.;
0294
0295   double m_mvtxDevs[6] = {0};
0296   double m_inttDevs[6] = {0};
0297   double m_tpcDevs[6] = {0};
0298   double m_mmDevs[6] = {0};
0299
0300   bool mvtxParam = false;
0301   bool inttParam = false;
0302   bool tpcParam = false;
0303   bool mmParam = false;
0304 };
0305
0306 #endif