packages/TrackerMillepedeAlignment/HelicalFitter.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef HELICALFITTER_H
0004 #define HELICALFITTER_H
0005
0006 #include "AlignmentDefs.h"
0007
0008 #include <tpc/TpcGlobalPositionWrapper.h>
0009
0010 #include <trackbase/ActsGeometry.h>
0011 #include <trackbase/ClusterErrorPara.h>
0012 #include <trackbase/TrackFitUtils.h>
0013
0014 #include <phparameter/PHParameterInterface.h>
0015 #include <tpc/TpcClusterZCrossingCorrection.h>
0016
0017 #include <fun4all/SubsysReco.h>
0018
0019 #include <map>
0020 #include <string>
0021
0022 class TpcClusterZCrossingCorrection;
0023 class PHCompositeNode;
0024 class TrackSeedContainer;
0025 class TrackSeed;
0026 class TrkrClusterContainer;
0027 class TF1;
0028 class TNtuple;
0029 class TFile;
0030 class Mille;
0031 class SvtxTrackSeed;
0032 class SvtxTrackMap;
0033 class SvtxVertexMap;
0034 class SvtxAlignmentStateMap;
0035 class SvtxTrack;
0036
0037 class HelicalFitter : public SubsysReco, public PHParameterInterface
0038 {
0039  public:
0040   HelicalFitter(const std::string& name = "HelicalFitter");
0041
0042   void SetDefaultParameters() override;
0043
0044   void set_field_dir(const double rescale)
0045   {
0046     _fieldDir = -1;
0047     if (rescale > 0)
0048     {
0049       _fieldDir = 1;
0050     }
0051   }
0052   void set_field(const std::string& field) { _field = field; }
0053
0054   int InitRun(PHCompositeNode* topNode) override;
0055
0056   int process_event(PHCompositeNode*) override;
0057
0058   int End(PHCompositeNode*) override;
0059
0060   void set_silicon_track_map_name(const std::string& map_name) { _silicon_track_map_name = map_name; }
0061   void set_track_map_name(const std::string& map_name) { _track_map_name = map_name; }
0062
0063   void set_use_event_vertex(bool flag) { use_event_vertex = flag; }
0064   void set_datafile_name(const std::string& file) { data_outfilename = file; }
0065   void set_steeringfile_name(const std::string& file) { steering_outfilename = file; }
0066   void set_mvtx_grouping(int group) { mvtx_grp = (AlignmentDefs::mvtxGrp) group; }
0067   void set_intt_grouping(int group) { intt_grp = (AlignmentDefs::inttGrp) group; }
0068   void set_tpc_grouping(int group) { tpc_grp = (AlignmentDefs::tpcGrp) group; }
0069   void set_mms_grouping(int group) { mms_grp = (AlignmentDefs::mmsGrp) group; }
0070   void set_test_output(bool test) { test_output = test; }
0071   void set_intt_layer_fixed(unsigned int layer);
0072   void set_mvtx_layer_fixed(unsigned int layer, unsigned int clamshell);
0073   void set_tpc_sector_fixed(unsigned int region, unsigned int sector, unsigned int side);
0074   void set_layer_param_fixed(unsigned int layer, unsigned int param);
0075   void set_ntuplefile_name(const std::string& file) { ntuple_outfilename = file; }
0076   void set_vertex_param_fixed(unsigned int param){ fixed_vertex_params.insert(param);}
0077   void set_straight_line_fit(bool flag) {straight_line_fit = flag; }
0078   void set_eta_cut(double eta_cut) {m_eta_cut = eta_cut;}
0079   //-1 is regular operation, 0 is east fixed, 1 is west fixed
0080   void set_do_mvtx_half(int half) {do_mvtx_half = half; }
0081   void set_fitted_subsystems(bool si, bool tpc, bool full)
0082   {
0083     fitsilicon = si;
0084     fittpc = tpc;
0085     fitfulltrack = full;
0086   }
0087
0088   void set_error_inflation_factor(unsigned int layer, float factor)
0089   {
0090     _layerMisalignment.insert(std::make_pair(layer, factor));
0091   }
0092
0093   void set_vertex(Acts::Vector3 inVertex, Acts::Vector3 inVertex_uncertainty)
0094   {
0095     vertexPosition = inVertex;
0096     vertexPosUncertainty = inVertex_uncertainty;
0097   }
0098
0099   void set_vtx_sigma(float x_sigma, float y_sigma)
0100   {
0101     vtx_sigma(0) = x_sigma;
0102     vtx_sigma(1) = y_sigma;
0103   }
0104
0105   // utility functions for analysis modules
0106   std::vector<float> fitClusters(std::vector<Acts::Vector3>& global_vec, std::vector<TrkrDefs::cluskey> cluskey_vec);
0107
0108   void getTrackletClusters(TrackSeed* _track, std::vector<Acts::Vector3>& global_vec, std::vector<TrkrDefs::cluskey>& cluskey_vec);
0109   Acts::Vector3 get_helix_pca(std::vector<float>& fitpars, const Acts::Vector3& global);
0110   void correctTpcGlobalPositions(std::vector<Acts::Vector3> global_vec, const std::vector<TrkrDefs::cluskey> &cluskey_vec);
0111   unsigned int addSiliconClusters(std::vector<float>& fitpars, std::vector<Acts::Vector3>& global_vec, std::vector<TrkrDefs::cluskey>& cluskey_vec);
0112
0113   void set_dca_cut(float dca) { dca_cut = dca; }
0114
0115  private:
0116   Mille* _mille;
0117
0118   int GetNodes(PHCompositeNode* topNode);
0119   int CreateNodes(PHCompositeNode* topNode);
0120   void getTrackletClusterList(TrackSeed* tracklet, std::vector<TrkrDefs::cluskey>& cluskey_vec);
0121
0122   Acts::Vector3 getPCALinePoint(const Acts::Vector3& global, const Acts::Vector3& tangent, const Acts::Vector3& posref);
0123   Acts::Vector3 get_line_plane_intersection(const Acts::Vector3& PCA, const Acts::Vector3& tangent,
0124                                             const Acts::Vector3& sensor_center, const Acts::Vector3& sensor_normal);
0125   std::pair<Acts::Vector3, Acts::Vector3> get_helix_tangent(const std::vector<float>& fitpars, Acts::Vector3 global);
0126   Acts::Vector3 get_helix_surface_intersection(const Surface& surf, std::vector<float>& fitpars, Acts::Vector3 global);
0127   Acts::Vector3 get_helix_surface_intersection(const Surface& surf, std::vector<float>& fitpars, Acts::Vector3 global, Acts::Vector3& pca, Acts::Vector3& tangent);
0128
0129   Acts::Vector3 get_helix_vtx(Acts::Vector3 event_vtx, const std::vector<float>& fitpars);
0130   Acts::Vector3 get_line_vtx(Acts::Vector3 event_vtx, const std::vector<float>& fitpars);
0131
0132   float convertTimeToZ(TrkrDefs::cluskey cluster_key, TrkrCluster* cluster);
0133   void makeTpcGlobalCorrections(TrkrDefs::cluskey cluster_key, short int crossing, Acts::Vector3& global);
0134
0135   Acts::Vector2 getClusterError(TrkrCluster* cluster, TrkrDefs::cluskey cluskey, Acts::Vector3& global);
0136
0137   bool is_tpc_sector_fixed(unsigned int layer, unsigned int sector, unsigned int side);
0138   bool is_mvtx_layer_fixed(unsigned int layer, unsigned int stave);
0139   bool is_intt_layer_fixed(unsigned int layer);
0140   bool is_layer_param_fixed(unsigned int layer, unsigned int param);
0141   bool is_vertex_param_fixed(unsigned int param);
0142
0143   void getLocalDerivativesXY(const Surface& surf, const Acts::Vector3& global, const std::vector<float>& fitpars, float lcl_derivativeX[5], float lcl_derivativeY[5], unsigned int layer);
0144   void getLocalDerivativesZeroFieldXY(const Surface& surf,  const Acts::Vector3& global, const std::vector<float>& fitpars, float lcl_derivativeX[5], float lcl_derivativeY[5], unsigned int layer);
0145
0146   void getLocalVtxDerivativesXY(SvtxTrack& track, const Acts::Vector3& track_vtx, const std::vector<float>& fitpars, float lcl_derivativeX[5], float lcl_derivativeY[5]);
0147   void getLocalVtxDerivativesZeroFieldXY(SvtxTrack& track, const Acts::Vector3& event_vtx, const std::vector<float>& fitpars, float lcl_derivativeX[5], float lcl_derivativeY[5]);
0148
0149   void getGlobalDerivativesXY(const Surface& surf, const Acts::Vector3& global, const Acts::Vector3& fitpoint, const std::vector<float>& fitpars, float glb_derivativeX[6], float glbl_derivativeY[6], unsigned int layer);
0150
0151   void getGlobalVtxDerivativesXY(SvtxTrack& track, const Acts::Vector3& track_vtx, float glbl_derivativeX[3], float glbl_derivativeY[3]);
0152
0153   void get_projectionXY(const Surface& surf, const std::pair<Acts::Vector3, Acts::Vector3>& tangent, Acts::Vector3& projX, Acts::Vector3& projY);
0154   void get_projectionVtxXY(SvtxTrack& track, const Acts::Vector3& event_vtx, Acts::Vector3& projX, Acts::Vector3& projY);
0155
0156   float getVertexResidual(Acts::Vector3 vtx);
0157
0158   void get_dca(SvtxTrack& track, float& dca3dxy, float& dca3dz, float& dca3dxysigma, float& dca3dzsigma, const Acts::Vector3& vertex);
0159   void get_dca_zero_field(SvtxTrack& track, float& dca3dxy, float& dca3dz, float& dca3dxysigma, float& dca3dzsigma, const Acts::Vector3& event_vertex);
0160
0161   std::pair<Acts::Vector3, Acts::Vector3> get_line(const std::vector<float>& fitpars);
0162   std::pair<Acts::Vector3, Acts::Vector3> get_line_zero_field(const std::vector<float>& fitpars);
0163   std::pair<Acts::Vector3, Acts::Vector3> get_line_tangent(const std::vector<float>& fitpars, Acts::Vector3 global);
0164   Acts::Vector3 get_line_surface_intersection(const Surface& surf, std::vector<float>& fitpars);
0165   Acts::Vector3 globalvtxToLocalvtx(SvtxTrack& track, const Acts::Vector3& event_vertex);
0166   Acts::Vector3 globalvtxToLocalvtx(SvtxTrack& track, const Acts::Vector3& event_vertex, Acts::Vector3 PCA);
0167   Acts::Vector3 localvtxToGlobalvtx(SvtxTrack& track, const Acts::Vector3& event_vtx, const Acts::Vector3& PCA);
0168
0169   //! global position wrapper
0170   TpcGlobalPositionWrapper m_globalPositionWrapper;
0171
0172   bool test_output = false;
0173
0174   ClusterErrorPara _ClusErrPara;
0175
0176   std::set<std::pair<unsigned int, unsigned int>> fixed_mvtx_layers;
0177   std::set<unsigned int> fixed_intt_layers;
0178   std::set<unsigned int> fixed_sectors;
0179   std::set<std::pair<unsigned int, unsigned int>> fixed_layer_params;
0180   std::set<unsigned int> fixed_vertex_params;
0181
0182   // set default groups to lowest level
0183   AlignmentDefs::mvtxGrp mvtx_grp = AlignmentDefs::mvtxGrp::snsr;
0184   AlignmentDefs::inttGrp intt_grp = AlignmentDefs::inttGrp::chp;
0185   AlignmentDefs::tpcGrp tpc_grp = AlignmentDefs::tpcGrp::htst;
0186   AlignmentDefs::mmsGrp mms_grp = AlignmentDefs::mmsGrp::tl;
0187
0188   //  TrackSeedContainer *_svtx_seed_map{nullptr};
0189   TrackSeedContainer* _track_map_tpc{nullptr};
0190   TrackSeedContainer* _track_map_silicon{nullptr};
0191   TrkrClusterContainer* _cluster_map{nullptr};
0192   ActsGeometry* _tGeometry{nullptr};
0193
0194   std::string data_outfilename{"mille_helical_output_data_file.bin"};
0195   std::string steering_outfilename{"steer_helical.txt"};
0196   std::string ntuple_outfilename{"HF_ntuple.root"};
0197
0198   TpcClusterZCrossingCorrection m_clusterCrossingCorrection;
0199
0200   bool fitsilicon{true};
0201   bool fittpc{false};
0202   bool fitfulltrack{false};
0203
0204   float dca_cut{0.19};  // cm
0205
0206   float m_eta_cut{99999.};
0207
0208   SvtxVertexMap* m_vertexmap{nullptr};
0209   SvtxTrackMap* m_trackmap{nullptr};
0210   SvtxAlignmentStateMap* m_alignmentmap{nullptr};
0211
0212   std::string _field;
0213   int _fieldDir{-1};
0214   std::map<unsigned int, float> _layerMisalignment;
0215
0216   std::string _track_map_name{"TpcTrackSeedContainer"};
0217   std::string _silicon_track_map_name{"SiliconTrackSeedContainer"};
0218
0219   bool make_ntuple{true};
0220   TNtuple* ntp{nullptr};
0221   TNtuple* track_ntp{nullptr};
0222   TFile* fout{nullptr};
0223
0224   bool use_event_vertex{false};
0225   bool use_intt_zfit{false};
0226   bool straight_line_fit = false;
0227   int do_mvtx_half = -1;
0228
0229   int event{0};
0230
0231   Acts::Vector3 vertexPosition;
0232   Acts::Vector3 vertexPosUncertainty;
0233   Acts::Vector2 vtx_sigma;
0234 };
0235
0236 #endif  // HELICALFITTER_H