packages/trackreco/PHSiliconTpcTrackMatching.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef PHSILICONTPCTRACKMATCHING_H
0004 #define PHSILICONTPCTRACKMATCHING_H
0005
0006 #include <fun4all/SubsysReco.h>
0007 #include <phparameter/PHParameterInterface.h>
0008 #include <tpc/TpcClusterZCrossingCorrection.h>
0009 #include <trackbase/ActsGeometry.h>
0010
0011 #include <map>
0012 #include <string>
0013
0014 class PHCompositeNode;
0015 class TrackSeedContainer;
0016 class TrackSeed;
0017 class TrkrClusterContainer;
0018 class TF1;
0019 class TrkrClusterCrossingAssoc;
0020 class TFile;
0021 class TNtuple;
0022
0023 class PHSiliconTpcTrackMatching : public SubsysReco, public PHParameterInterface
0024 {
0025  public:
0026   PHSiliconTpcTrackMatching(const std::string &name = "PHSiliconTpcTrackMatching");
0027
0028   ~PHSiliconTpcTrackMatching() override;
0029
0030   void SetDefaultParameters() override;
0031
0032   void set_phi_search_window(const double win) { _phi_search_win = win; }
0033   void set_eta_search_window(const double win) { _eta_search_win = win; }
0034   void set_x_search_window(const double win) { _x_search_win = win; }
0035   void set_y_search_window(const double win) { _y_search_win = win; }
0036   void set_z_search_window(const double win) { _z_search_win = win; }
0037   void set_crossing_deltaz_max(const double dz) {_crossing_deltaz_max = dz;}
0038   void set_crossing_deltaz_min(const double dz) {_crossing_deltaz_min = dz;}
0039   void set_deltaeta_min(const double deta) {_deltaeta_min = deta;}
0040
0041   float get_phi_search_window() const { return _phi_search_win; }
0042   float get_eta_search_window() const { return _eta_search_win; }
0043   float get_x_search_window() const { return _x_search_win; }
0044   float get_y_search_window() const { return _y_search_win; }
0045   float get_z_search_window() const { return _z_search_win; }
0046
0047   // 2024/01/22 update
0048   struct WindowMatcher {
0049     // --- new method, comparing to a+b*exp(c/pT)
0050     // Each Arr3D object contains a,b,c in order
0051     // Up to four curved are needed:
0052     // - for for positive and negative tracks
0053     // - only one curve if |dX|<fn_max, or two if like  fn_min<dX<fnmax
0054     using Arr3D = std::array<double,3>;
0055     std::string print_fn(const Arr3D&);
0056     Arr3D posLo { 100, 0, 0 }; // (above a2,b2,c2), 100 for |dX|
0057     Arr3D posHi { 100, 0, 0 }; // (above a3,b3,c3), 100 for use _use_legacy_windowing
0058     Arr3D negLo { 100, 0, 0 }; // (above a0,b0,c0), 100 for |dX|
0059     Arr3D negHi { 100, 0, 0 }; // (above a1,b1,c1), 100 for treat all tracks pos Q
0060
0061     // efficiency flags set during PHSiliconTpcTrackMatching::InitRun()
0062     bool fabs_max_posQ  = true;
0063     bool fabs_max_negQ  = true;
0064     bool negLo_b0 = true;
0065     bool negHi_b0 = true;
0066     bool posLo_b0 = true;
0067     bool posHi_b0 = true;
0068     double min_pt_posQ  = 0.25; // only grow function windows down to 150 MeV
0069     double min_pt_negQ  = 0.25; // only grow function windows down to 150 MeV
0070
0071     WindowMatcher(
0072         const Arr3D& _posLo={100,0,0},
0073         const Arr3D& _posHi={100,0,0},
0074         const Arr3D& _negLo={100,0,0},
0075         const Arr3D& _negHi={100,0,0},
0076         const double _min_pt_posQ=0.25,
0077         const double _min_pt_negQ=0.25)
0078       : posLo{_posLo}, posHi{_posHi}, negLo{_negLo}, negHi{_negHi},
0079       min_pt_posQ{_min_pt_posQ}, min_pt_negQ{_min_pt_negQ} {};
0080
0081     inline double fn_exp(const Arr3D& arr, const bool& b_is_0, double pT) {
0082       return (b_is_0 ? arr[0] : arr[0]+arr[1]*exp(arr[2]/pT));
0083     }
0084
0085     void init_bools(const std::string& which_window="", const bool print=false);
0086
0087     bool in_window(bool posQ, const double tpc_pt, const double tpc_X, const double si_X);
0088
0089     // initialize to fn_lo < deltaX < fn_hi for +Q, and fn_lo < deltaX < fn_hi for -Q
0090
0091     void reset_fns() {
0092       posLo={100,0,0};
0093       posHi={100,0,0};
0094       negLo={100,0,0};
0095       negHi={100,0,0};
0096     };
0097
0098     // same max for |deltaX| for pos and neg Q
0099     void set_QoverpT_maxabs    (const Arr3D& _posHi, const double _min_pt=0.25)
0100     { reset_fns(); posHi=_posHi; min_pt_posQ = _min_pt; };
0101
0102     // same range for deltaX for pos and neg Q
0103     void set_QoverpT_range     (const Arr3D& _posLo, const Arr3D& _posHi, const double _min_pt=0.25)
0104     { reset_fns(); posLo=_posLo; posHi=_posHi; min_pt_posQ=_min_pt; };
0105
0106     // max for |deltaX| for pos Q
0107     void set_posQoverpT_maxabs (const Arr3D& _posHi, const double _min_pt=0.25)
0108     { posLo={100.,0.,0.}; posHi=_posHi; min_pt_posQ = _min_pt; };
0109
0110     // max for |deltaX| for neg Q
0111     void set_negQoverpT_maxabs (const Arr3D& _negHi, const double _min_pt=0.25)
0112     { posLo={100.,0.,0.}; negHi=_negHi; min_pt_negQ = _min_pt; };
0113
0114     // range for deltaX for pos Q
0115     void set_posQoverpT_range  (const Arr3D& _posLo, const Arr3D& _posHi, const double _min_pt=0.25)
0116     { posLo=_posLo; posHi=_posHi; min_pt_posQ = _min_pt; };
0117
0118     // range for deltaX for neg Q
0119     void set_negQoverpT_range  (const Arr3D& _negLo, const Arr3D& _negHi, const double _min_pt=0.25)
0120     { negLo=_negLo; negHi=_negHi; min_pt_negQ = _min_pt; };
0121
0122   };
0123
0124   // initialize the window matchers with default values
0125   WindowMatcher window_dx   { {100.,0.,0.}, {5.3, 0., 0.} };
0126   WindowMatcher window_dy   { {100.,0.,0.}, {5.2, 0., 0.} };
0127   WindowMatcher window_dz   { {100.,0.,0.}, {0., 2.6, 0.38}, {100,0.,0.,}, {0., 1.45, 0.49} };
0128   WindowMatcher window_dphi { {-0.25, 0., 0.},  {0.05, 0., 0.} };
0129   WindowMatcher window_deta { {100.,0.,0.}, {0.050, 0.0064, 1.1}, {100,0.,0.,}, {0.045, 0.0031, 1.0} };
0130
0131   bool _print_windows = false;
0132   void print_windows(bool print=true) { _print_windows = print; }
0133
0134   void zeroField(const bool flag) { _zero_field = flag; }
0135
0136   //  void set_use_old_matching(const bool flag) { _use_old_matching = flag; }
0137
0138   void set_test_windows_printout(const bool test) { _test_windows = test; }
0139   void set_file_name(const std::string &name) { _file_name = name; }
0140   void set_pp_mode(const bool flag) { _pp_mode = flag; }
0141   void set_use_intt_crossing(const bool flag) { _use_intt_crossing = flag; }
0142
0143   int InitRun(PHCompositeNode *topNode) override;
0144
0145   int process_event(PHCompositeNode *) override;
0146
0147   int End(PHCompositeNode *) override;
0148
0149   void fieldMap(std::string &fieldmap) { m_fieldMap = fieldmap; }
0150
0151   void set_silicon_track_map_name(const std::string &map_name) { _silicon_track_map_name = map_name; }
0152   void set_track_map_name(const std::string &map_name) { _track_map_name = map_name; }
0153   void SetIteration(int iter) { _n_iteration = iter; }
0154
0155  private:
0156   int GetNodes(PHCompositeNode *topNode);
0157
0158   void findEtaPhiMatches(std::set<unsigned int> &tpc_matched_set,
0159                          std::set<unsigned int> &tpc_unmatched_set,
0160                          std::multimap<unsigned int, unsigned int> &tpc_matches);
0161   std::vector<short int> getInttCrossings(TrackSeed *si_track);
0162   void checkZMatches(std::multimap<unsigned int, unsigned int> &tpc_matches,
0163              std::multimap<unsigned int, unsigned int> &bad_map);
0164   short int getCrossingIntt(TrackSeed *_tracklet_si);
0165   // void findCrossingGeometrically(std::multimap<unsigned int, unsigned int> tpc_matches);
0166   short int findCrossingGeometrically(unsigned int tpc_id, unsigned int si_id);
0167   double getBunchCrossing(unsigned int trid, double z_mismatch);
0168
0169   TFile *_file = nullptr;
0170   TNtuple *_tree = nullptr;
0171
0172   std::string _file_name = "track_match.root";
0173
0174   // default values, can be replaced from the macro
0175   double _phi_search_win = 0.01;
0176   double _eta_search_win = 0.004;
0177   double _x_search_win = 0.3;
0178   double _y_search_win = 0.3;
0179   double _z_search_win = 0.4;
0180
0181   //  bool _use_old_matching = false;  // normally false
0182
0183   bool _zero_field = false;     // fit straight lines if true
0184
0185   TrackSeedContainer *_svtx_seed_map{nullptr};
0186   TrackSeedContainer *_track_map{nullptr};
0187   TrackSeedContainer *_track_map_silicon{nullptr};
0188   TrackSeed *_tracklet_tpc{nullptr};
0189   TrackSeed *_tracklet_si{nullptr};
0190   TrkrClusterContainer *_cluster_map{nullptr};
0191   ActsGeometry *_tGeometry{nullptr};
0192   TrkrClusterCrossingAssoc *_cluster_crossing_map{nullptr};
0193   int m_event = 0;
0194   std::map<unsigned int, double> _z_mismatch_map;
0195
0196   TpcClusterZCrossingCorrection _clusterCrossingCorrection;
0197   float _crossing_deltaz_max = 10.0;
0198   float _crossing_deltaz_min = 1.5;
0199   float _deltaeta_min = 0.03;
0200
0201   //  double _collision_rate = 50e3;  // input rate for phi correction
0202   //  double _reference_collision_rate = 50e3;  // reference rate for phi correction
0203   //  double _si_vertex_dzmax = 0.25;  // mm
0204   double fieldstrength{std::numeric_limits<double>::quiet_NaN()};
0205
0206   bool _test_windows = false;
0207   bool _pp_mode = false;
0208   bool _use_intt_crossing = true;  // should always be true except for testing
0209
0210   int _n_iteration = 0;
0211   std::string _track_map_name = "TpcTrackSeedContainer";
0212   std::string _silicon_track_map_name = "SiliconTrackSeedContainer";
0213   std::string m_fieldMap = "1.4";
0214   std::vector<TrkrDefs::cluskey> getTrackletClusterList(TrackSeed* tracklet);
0215 };
0216
0217 #endif  //  PHSILICONTPCTRACKMATCHING_H