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 PHG4TpcGeomContainer;
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
0133   void set_svtx_alignment_state_map_name(const std::string& map_name) {
0134       _svtx_alignment_state_map_name = map_name;
0135       m_alignStates.alignmentStateMap(map_name);
0136   }
0137
0138   /// Set flag for pp running
0139   void set_pp_mode(bool ispp) { m_pp_mode = ispp; }
0140
0141   void set_enable_geometric_crossing_estimate(bool flag) { m_enable_crossing_estimate = flag ; }
0142   void set_use_clustermover(bool use) { m_use_clustermover = use; }
0143   void ignoreLayer(int layer) { m_ignoreLayer.insert(layer); }
0144   void setTrkrClusterContainerName(std::string &name){ m_clusterContainerName = name; }
0145   void setDirectNavigation(bool flag) { m_directNavigation = flag; }
0146
0147  private:
0148   /// Get all the nodes
0149   int getNodes(PHCompositeNode* topNode);
0150
0151   /// Create new nodes
0152   int createNodes(PHCompositeNode* topNode);
0153
0154   void loopTracks(Acts::Logging::Level logLevel);
0155
0156   /// Convert the acts track fit result to an svtx track
0157   void updateSvtxTrack(
0158     const std::vector<Acts::MultiTrajectoryTraits::IndexType>& tips,
0159     const Trajectory::IndexedParameters& paramsMap,
0160     const 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(const FitResult& fitOutput, TrackSeed* seed,
0181                          SvtxTrack* track,
0182                          const 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   //! tracks
0256 //  SvtxTrackMap* m_seedTracks = nullptr;
0257
0258   //! tpc global position wrapper
0259   TpcGlobalPositionWrapper m_globalPositionWrapper;
0260
0261   //! list of layers to be removed from fit
0262   std::set<int> m_ignoreLayer;
0263
0264   bool m_use_clustermover = true;
0265
0266   std::string m_fieldMap;
0267
0268   int _n_iteration = 0;
0269   std::string _track_map_name = "SvtxTrackMap";
0270   std::string _svtx_seed_map_name = "SvtxTrackSeedContainer";
0271   std::string _svtx_alignment_state_map_name =  "SvtxAlignmentStateMap";
0272
0273   /// Default particle assumption to pion
0274   unsigned int m_pHypothesis = 211;
0275
0276   SvtxAlignmentStateMap* m_alignmentStateMap = nullptr;
0277   ActsAlignmentStates m_alignStates;
0278   bool m_commissioning = false;
0279
0280   bool m_enable_crossing_estimate = false;
0281
0282   PHG4TpcGeomContainer* _tpccellgeo = nullptr;
0283
0284   /// Variables for doing event time execution analysis
0285   bool m_timeAnalysis = false;
0286   TFile* m_timeFile = nullptr;
0287   TH1* h_eventTime = nullptr;
0288   TH2* h_fitTime = nullptr;
0289   TH1* h_updateTime = nullptr;
0290   TH1* h_stateTime = nullptr;
0291   TH1* h_rotTime = nullptr;
0292
0293   std::vector<const Acts::Surface*> m_materialSurfaces = {};
0294
0295   struct MaterialSurfaceSelector {
0296     std::vector<const Acts::Surface*> surfaces = {};
0297
0298     /// @param surface is the test surface
0299     void operator()(const Acts::Surface* surface) {
0300       if (surface->surfaceMaterial() != nullptr) {
0301         if (std::find(surfaces.begin(), surfaces.end(), surface) ==
0302             surfaces.end()) {
0303           surfaces.push_back(surface);
0304         }
0305       }
0306     }
0307   };
0308 };
0309
0310 #endif