packages/trackreco/PHActsTrkFitter.h

0001 /*!
0002  *  \file       PHActsTrkFitter.h
0003  *  \brief      Refit SvtxTracks with Acts.
0004  *  \details            Refit SvtxTracks with Acts
0005  *  \author     Joe Osborn, Tony Frawley <afrawley@fsu.edu>
0006  */
0007
0008 #ifndef TRACKRECO_ACTSTRKFITTER_H
0009 #define TRACKRECO_ACTSTRKFITTER_H
0010
0011 #include "ActsAlignmentStates.h"
0012 #include "ActsEvaluator.h"
0013
0014 #include <fun4all/SubsysReco.h>
0015
0016 #include <trackbase/ActsSourceLink.h>
0017 #include <trackbase/ActsTrackFittingAlgorithm.h>
0018
0019 #include <tpc/TpcGlobalPositionWrapper.h>
0020
0021 #include <Acts/Definitions/Algebra.hpp>
0022 #include <Acts/EventData/VectorMultiTrajectory.hpp>
0023 #include <Acts/Utilities/BinnedArray.hpp>
0024 #include <Acts/Utilities/Logger.hpp>
0025 #include <Acts/Utilities/Helpers.hpp>
0026
0027 #include <ActsExamples/EventData/Trajectories.hpp>
0028
0029 #include <TFile.h>
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <memory>
0033 #include <string>
0034
0035 class alignmentTransformationContainer;
0036 class ActsGeometry;
0037 class SvtxTrack;
0038 class SvtxTrackMap;
0039 class TrackSeed;
0040 class TrackSeedContainer;
0041 class TrkrClusterContainer;
0042 class SvtxAlignmentStateMap;
0043 class PHG4TpcCylinderGeomContainer;
0044
0045 using SourceLink = ActsSourceLink;
0046 using FitResult = ActsTrackFittingAlgorithm::TrackFitterResult;
0047 using Trajectory = ActsExamples::Trajectories;
0048 using Measurement = Acts::Measurement<Acts::BoundIndices, 2>;
0049 using SurfacePtrVec = std::vector<const Acts::Surface*>;
0050 using SourceLinkVec = std::vector<Acts::SourceLink>;
0051
0052 class PHActsTrkFitter : public SubsysReco
0053 {
0054  public:
0055   /// Default constructor
0056   PHActsTrkFitter(const std::string& name = "PHActsTrkFitter");
0057
0058   /// Destructor
0059   ~PHActsTrkFitter() override = default;
0060
0061   /// End, write and close files
0062   int End(PHCompositeNode* topNode) override;
0063
0064   /// Get and create nodes
0065   int InitRun(PHCompositeNode* topNode) override;
0066
0067   /// Process each event by calling the fitter
0068   int process_event(PHCompositeNode* topNode) override;
0069
0070   int ResetEvent(PHCompositeNode* topNode) override;
0071
0072   /// Do some internal time benchmarking analysis
0073   void doTimeAnalysis(bool timeAnalysis) { m_timeAnalysis = timeAnalysis; }
0074
0075   /// Run the direct navigator to fit only tracks with silicon+MM hits
0076   void fitSiliconMMs(bool fitSiliconMMs)
0077   {
0078     m_fitSiliconMMs = fitSiliconMMs;
0079   }
0080
0081   /// with direct navigation, force a fit with only silicon hits
0082   void forceSiOnlyFit(bool forceSiOnlyFit)
0083   {
0084     m_forceSiOnlyFit = forceSiOnlyFit;
0085   }
0086
0087   /// require micromegas in SiliconMM fits
0088   void setUseMicromegas(bool value)
0089   {
0090     m_useMicromegas = value;
0091   }
0092   void ignoreSilicon()
0093   {
0094     m_ignoreSilicon = true;
0095   }
0096   void setUpdateSvtxTrackStates(bool fillSvtxTrackStates)
0097   {
0098     m_fillSvtxTrackStates = fillSvtxTrackStates;
0099   }
0100
0101   void useActsEvaluator(bool actsEvaluator)
0102   {
0103     m_actsEvaluator = actsEvaluator;
0104   }
0105   void useActsEvaluatorSimulation(bool actsEvaluator)
0106   {
0107     m_simActsEvaluator = actsEvaluator;
0108     m_actsEvaluator = actsEvaluator;
0109   }
0110   void setEvaluatorName(const std::string& name) { m_evalname = name; }
0111   void setFieldMap(const std::string& fieldMap)
0112   {
0113     m_fieldMap = fieldMap;
0114   }
0115
0116   void setAbsPdgHypothesis(unsigned int pHypothesis)
0117   {
0118     m_pHypothesis = pHypothesis;
0119   }
0120
0121   void commissioning(bool com) { m_commissioning = com; }
0122
0123   void useOutlierFinder(bool outlier) { m_useOutlierFinder = outlier; }
0124   void setOutlierFinderOutfile(const std::string& outfilename)
0125   {
0126     m_outlierFinder.outfileName(outfilename);
0127   }
0128
0129   void SetIteration(int iter) { _n_iteration = iter; }
0130   void set_track_map_name(const std::string& map_name) { _track_map_name = map_name; }
0131   void set_svtx_seed_map_name(const std::string& map_name) { _svtx_seed_map_name = map_name; }
0132   void set_trajctories_name(const std::string& map_name) {m_trajectories_name = map_name; }
0133
0134   void set_svtx_alignment_state_map_name(const std::string& map_name) {
0135       _svtx_alignment_state_map_name = map_name;
0136       m_alignStates.alignmentStateMap(map_name);
0137   }
0138
0139   /// Set flag for pp running
0140   void set_pp_mode(bool ispp) { m_pp_mode = ispp; }
0141
0142   void set_enable_geometric_crossing_estimate(bool flag) { m_enable_crossing_estimate = flag ; }
0143   void set_use_clustermover(bool use) { m_use_clustermover = use; }
0144   void ignoreLayer(int layer) { m_ignoreLayer.insert(layer); }
0145   void setTrkrClusterContainerName(std::string &name){ m_clusterContainerName = name; }
0146   void setDirectNavigation(bool flag) { m_directNavigation = flag; }
0147
0148  private:
0149   /// Get all the nodes
0150   int getNodes(PHCompositeNode* topNode);
0151
0152   /// Create new nodes
0153   int createNodes(PHCompositeNode* topNode);
0154
0155   void loopTracks(Acts::Logging::Level logLevel);
0156
0157   /// Convert the acts track fit result to an svtx track
0158   void updateSvtxTrack(std::vector<Acts::MultiTrajectoryTraits::IndexType>& tips,
0159                        Trajectory::IndexedParameters& paramsMap,
0160                        ActsTrackFittingAlgorithm::TrackContainer& tracks,
0161                        SvtxTrack* track);
0162
0163   /// Helper function to call either the regular navigation or direct
0164   /// navigation, depending on m_fitSiliconMMs
0165   ActsTrackFittingAlgorithm::TrackFitterResult fitTrack(
0166       const std::vector<Acts::SourceLink>& sourceLinks,
0167       const ActsTrackFittingAlgorithm::TrackParameters& seed,
0168       const ActsTrackFittingAlgorithm::GeneralFitterOptions&
0169           kfOptions,
0170       const SurfacePtrVec& surfSequence,
0171       const CalibratorAdapter& calibrator,
0172       ActsTrackFittingAlgorithm::TrackContainer& tracks);
0173
0174   /// Functions to get list of sorted surfaces for direct navigation, if
0175   /// applicable
0176   SourceLinkVec getSurfaceVector(const SourceLinkVec& sourceLinks,
0177                                  SurfacePtrVec& surfaces) const;
0178   void checkSurfaceVec(SurfacePtrVec& surfaces) const;
0179
0180   bool getTrackFitResult(FitResult& fitOutput, TrackSeed* seed,
0181                          SvtxTrack* track,
0182                          ActsTrackFittingAlgorithm::TrackContainer& tracks,
0183                          const ActsTrackFittingAlgorithm::MeasurementContainer& measurements);
0184
0185   Acts::BoundSquareMatrix setDefaultCovariance() const;
0186   void printTrackSeed(const ActsTrackFittingAlgorithm::TrackParameters& seed) const;
0187
0188   /// Event counter
0189   int m_event = 0;
0190
0191   /// Options that Acts::Fitter needs to run from MakeActsGeometry
0192   ActsGeometry* m_tGeometry = nullptr;
0193
0194   /// Configuration containing the fitting function instance
0195   ActsTrackFittingAlgorithm::Config m_fitCfg;
0196
0197   /// TrackMap containing SvtxTracks
0198   alignmentTransformationContainer* m_alignmentTransformationMap = nullptr;  // added for testing purposes
0199   alignmentTransformationContainer* m_alignmentTransformationMapTransient = nullptr;
0200   std::set<Acts::GeometryIdentifier> m_transient_id_set;
0201   Acts::GeometryContext m_transient_geocontext;
0202   SvtxTrackMap* m_trackMap = nullptr;
0203   SvtxTrackMap* m_directedTrackMap = nullptr;
0204   TrkrClusterContainer* m_clusterContainer = nullptr;
0205   TrackSeedContainer* m_seedMap = nullptr;
0206   TrackSeedContainer* m_tpcSeeds = nullptr;
0207   TrackSeedContainer* m_siliconSeeds = nullptr;
0208
0209   /// Number of acts fits that returned an error
0210   int m_nBadFits = 0;
0211
0212   /// Boolean to use normal tracking geometry navigator or the
0213   /// Acts::DirectedNavigator with a list of sorted silicon+MM surfaces
0214   bool m_fitSiliconMMs = false;
0215
0216   bool m_forceSiOnlyFit = false;
0217
0218   /// requires micromegas present when fitting silicon-MM surfaces
0219   bool m_useMicromegas = true;
0220
0221   /// A bool to update the SvtxTrackState information (or not)
0222   bool m_fillSvtxTrackStates = true;
0223
0224   /// bool to ignore the silicon clusters in the fit
0225   bool m_ignoreSilicon = false;
0226
0227   /// A bool to use the chi2 outlier finder in the track fitting
0228   bool m_useOutlierFinder = false;
0229   ResidualOutlierFinder m_outlierFinder;
0230
0231   /// Flag for pp running
0232   bool m_pp_mode = false;
0233
0234   bool m_directNavigation = true;
0235
0236   // do we have a constant field
0237   bool m_ConstField{false};
0238   double fieldstrength{std::numeric_limits<double>::quiet_NaN()};
0239
0240   // max variation of bunch crossing away from crossing_estimate
0241   short int max_bunch_search = 2;
0242
0243   //name of TRKR_CLUSTER container
0244   std::string m_clusterContainerName = "TRKR_CLUSTER";
0245
0246   //!@name evaluator
0247   //@{
0248   bool m_actsEvaluator = false;
0249   bool m_simActsEvaluator = false;
0250   std::unique_ptr<ActsEvaluator>
0251       m_evaluator = nullptr;
0252   std::string m_evalname = "ActsEvaluator.root";
0253   //@}
0254
0255   //! acts trajectories
0256   std::map<const unsigned int, Trajectory>* m_trajectories = nullptr;
0257   std::string m_trajectories_name = "ActsTrajectories";
0258
0259   //! tracks
0260 //  SvtxTrackMap* m_seedTracks = nullptr;
0261
0262   //! tpc global position wrapper
0263   TpcGlobalPositionWrapper m_globalPositionWrapper;
0264
0265   //! list of layers to be removed from fit
0266   std::set<int> m_ignoreLayer;
0267
0268   bool m_use_clustermover = true;
0269
0270   std::string m_fieldMap;
0271
0272   int _n_iteration = 0;
0273   std::string _track_map_name = "SvtxTrackMap";
0274   std::string _svtx_seed_map_name = "SvtxTrackSeedContainer";
0275   std::string _svtx_alignment_state_map_name =  "SvtxAlignmentStateMap";
0276
0277   /// Default particle assumption to pion
0278   unsigned int m_pHypothesis = 211;
0279
0280   SvtxAlignmentStateMap* m_alignmentStateMap = nullptr;
0281   ActsAlignmentStates m_alignStates;
0282   bool m_commissioning = false;
0283
0284   bool m_enable_crossing_estimate = false;
0285
0286   PHG4TpcCylinderGeomContainer* _tpccellgeo = nullptr;
0287
0288   /// Variables for doing event time execution analysis
0289   bool m_timeAnalysis = false;
0290   TFile* m_timeFile = nullptr;
0291   TH1* h_eventTime = nullptr;
0292   TH2* h_fitTime = nullptr;
0293   TH1* h_updateTime = nullptr;
0294   TH1* h_stateTime = nullptr;
0295   TH1* h_rotTime = nullptr;
0296
0297   std::vector<const Acts::Surface*> m_materialSurfaces = {};
0298
0299   struct MaterialSurfaceSelector {
0300     std::vector<const Acts::Surface*> surfaces = {};
0301
0302     /// @param surface is the test surface
0303     void operator()(const Acts::Surface* surface) {
0304       if (surface->surfaceMaterial() != nullptr) {
0305         if (std::find(surfaces.begin(), surfaces.end(), surface) ==
0306             surfaces.end()) {
0307           surfaces.push_back(surface);
0308         }
0309       }
0310     }
0311   };
0312 };
0313
0314 #endif