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 PHG4TpcGeomContainer;
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   void set_sampa_tzero_bias(double tzb) { m_sampa_tzero_bias = tzb; }
0143   void set_apply_tpc_tzero_correction(bool flag) { m_apply_tpc_tzero_correction = flag; }
0144
0145   void set_nSurfPhi(unsigned int value)
0146   {
0147     m_nSurfPhi = value;
0148   }
0149   //  void set_maxSurfZ(double value) {m_maxSurfZ = value;}  // set to TPC gas volume length
0150
0151   void set_mvtx_applymisalign(bool b) { m_mvtxapplymisalign = b; }
0152   void set_intt_survey(bool surv) { m_inttSurvey = surv; }
0153
0154  private:
0155   /// Main function to build all acts geometry for use in the fitting modules
0156   int buildAllGeometry(PHCompositeNode *topNode);
0157
0158   //! Get all the nodes
0159   int getNodes(PHCompositeNode *);
0160
0161   //! Create New nodes
0162   int createNodes(PHCompositeNode *);
0163
0164   /// Functions to edit TGeoManager to include TPC boxes
0165   void setPlanarSurfaceDivisions();
0166   void editTPCGeometry(PHCompositeNode *topNode);
0167   void addActsTpcSurfaces(TGeoVolume *tpc_gas_vol,
0168                           TGeoManager *geoManager);
0169
0170   /// Silicon layers made by BuildSiliconLayers and its helper functions
0171   void buildActsSurfaces();
0172
0173   /// Function that mimics ActsExamples::GeometryExampleBase
0174   void makeGeometry(int argc, char *argv[],
0175                     ActsExamples::TGeoDetectorWithOptions &detector);
0176 #ifndef __CLING__
0177   std::pair<std::shared_ptr<const Acts::TrackingGeometry>,
0178             std::vector<std::shared_ptr<ActsExamples::IContextDecorator>>>
0179   build(const boost::program_options::variables_map &vm,
0180         ActsExamples::TGeoDetectorWithOptions &detector);
0181 #endif
0182   void readTGeoLayerBuilderConfigsFile(const std::string &path,
0183                                        ActsExamples::TGeoDetector::Config &config);
0184
0185   void setMaterialResponseFile(std::string &responseFile,
0186                                std::string &materialFile);
0187
0188   /// Get hitsetkey from TGeoNode for each detector geometry
0189   void getInttKeyFromNode(TGeoNode *gnode);
0190   void getMvtxKeyFromNode(TGeoNode *gnode);
0191   void getTpcKeyFromNode(TGeoNode *gnode);
0192
0193   /// Make the Surface<-->TrkrDef::hitsetkey map pairs for each of
0194   /// the various subdetectors
0195   void makeMvtxMapPairs(TrackingVolumePtr &mvtxVolume);
0196   void makeInttMapPairs(TrackingVolumePtr &inttVolume);
0197   void makeTpcMapPairs(TrackingVolumePtr &tpcVolume);
0198
0199   /// bind micromegas surfaces to hitset id
0200   void makeMmMapPairs(TrackingVolumePtr &tpcVolume);
0201
0202   /// Get subdetector hitsetkey from the local sensor unit coordinates
0203   TrkrDefs::hitsetkey getMvtxHitSetKeyFromCoords(unsigned int layer,
0204                                                  std::vector<double> &world);
0205   TrkrDefs::hitsetkey getInttHitSetKeyFromCoords(unsigned int layer,
0206                                                  std::vector<double> &world);
0207   TrkrDefs::hitsetkey getTpcHitSetKeyFromCoords(std::vector<double> &world);
0208
0209   //   /// Helper diagnostic function for identifying active layers in subdetectors
0210   //   void isActive(TGeoNode *gnode, int nmax_print);
0211
0212   //   /// Makes map of TrkrHitSetKey<-->TGeoNode
0213   //   void makeTGeoNodeMap(PHCompositeNode *topNode);
0214
0215   void unpackVolumes();
0216
0217   /// Subdetector geometry containers for getting layer information
0218   PHG4CylinderGeomContainer *m_geomContainerMvtx = nullptr;
0219   PHG4CylinderGeomContainer *m_geomContainerIntt = nullptr;
0220   PHG4CylinderGeomContainer *m_geomContainerMicromegas = nullptr;
0221   PHG4TpcGeomContainer *m_geomContainerTpc = nullptr;
0222   TGeoManager *m_geoManager = nullptr;
0223
0224   // Switch to use or not use the INTT survey geometry
0225   bool m_inttSurvey = true;
0226   const float m_inttbarrelcenter_survey_x = 0.4026857142857132 / 10.;
0227   const float m_inttbarrelcenter_survey_y = -2.886627321428573 / 10.;
0228
0229   // Switch for applying misalignment for mvtx
0230   bool m_mvtxapplymisalign = false;
0231   std::vector<double> v_globaldisplacement = {0., 0., 0.};
0232
0233   bool m_useField = true;
0234   std::map<uint8_t, double> m_misalignmentFactor;
0235
0236   /// Several maps that connect Acts world to sPHENIX G4 world
0237   std::map<TrkrDefs::hitsetkey, TGeoNode *> m_clusterNodeMap;
0238   std::map<TrkrDefs::hitsetkey, Surface> m_clusterSurfaceMapSilicon;
0239   std::map<unsigned int, std::vector<Surface>> m_clusterSurfaceMapTpcEdit;  // uses layer as key
0240   std::map<TrkrDefs::hitsetkey, Surface> m_clusterSurfaceMapMmEdit;
0241
0242   /// These don't change, we are building the tpc this way!
0243   static constexpr unsigned int m_nTpcLayers = 48;
0244   static constexpr unsigned int m_nTpcModulesPerLayer = 12;
0245   static constexpr unsigned int m_nTpcSides = 2;
0246
0247   /// TPC Acts::Surface subdivisions
0248   double m_minSurfZ = 0.;
0249   /// This value should be slightly less than the TPC gas volume in TGeo
0250   double m_maxSurfZ =  0;
0251   unsigned int m_nSurfZ = 1;
0252   unsigned int m_nSurfPhi = 12;
0253   double m_surfStepPhi = 0;
0254   double m_surfStepZ = 0;
0255   double m_moduleStepPhi = 0;
0256   double m_modulePhiStart = 0;
0257
0258   /// Debugger for printing out tpc active volumes
0259   //  int nprint_tpc = 0;
0260
0261   /// TPC TGeoManager editing box surfaces subdivisions
0262   const static int m_nTpcSectors = 3;
0263
0264   double m_layerRadius[m_nTpcLayers] = {0};
0265   double m_layerThickness[m_nTpcLayers] = {0};
0266
0267   // Spaces to prevent boxes from touching when placed
0268   const double half_width_clearance_thick = 0.4999;
0269   const double half_width_clearance_phi = 0.4999;
0270   /// z does not need spacing as the boxes are rotated around the z axis
0271   const double half_width_clearance_z = 0.5;
0272
0273   /// The acts geometry object
0274   ActsExamples::TGeoDetectorWithOptions m_detector;
0275
0276   /// Acts geometry objects that are needed to create (for example) the fitter
0277   TrackingGeometry m_tGeometry;
0278   std::shared_ptr<Acts::MagneticFieldProvider> m_magneticField;
0279   Acts::GeometryContext m_geoCtxt;
0280
0281   /// Structs to put on the node tree which carry around ActsGeom info
0282   ActsGeometry *m_actsGeometry = nullptr;
0283
0284   std::map<unsigned int, unsigned int> base_layer_map = {{10, 0}, {12, 3}, {14, 7}, {16, 55}};
0285   unsigned int mvtx_chips_per_stave = 9;
0286
0287   /// Verbosity value handed from PHActsSourceLinks
0288   //  int m_verbosity = 0;
0289
0290   double m_drift_velocity = 0.;  // cm/ns, override from macro
0291   double m_max_driftlength = 0.;  // override from macro
0292   double m_CM_halfwidth = 0.;  // central membrane half width in cm
0293
0294   bool m_apply_tpc_tzero_correction = false;
0295   double m_tpc_tzero = 0.0;  // ns, override from macro
0296   double m_sampa_tzero_bias = 0.0;  // ns, override from macro
0297
0298   /// Magnetic field components to set Acts magnetic field
0299   std::string m_magField = "1.4";
0300   double m_magFieldRescale = -1.;
0301
0302   double m_mvtxDevs[6] = {0};
0303   double m_inttDevs[6] = {0};
0304   double m_tpcDevs[6] = {0};
0305   double m_mmDevs[6] = {0};
0306
0307   bool mvtxParam = false;
0308   bool inttParam = false;
0309   bool tpcParam = false;
0310   bool mmParam = false;
0311 };
0312
0313 #endif