packages/trackreco/PHStreakFinder.h

0001 /*!
0002  *  \file RPHStreakFinder.h
0003  *  \RTree based hough tracking for cosmics
0004  *  \author Christof Roland
0005  */
0006
0007
0008 //begin
0009
0010 #include <fun4all/SubsysReco.h>
0011 #include <trackbase/TrkrDefs.h>  // for cluskey
0012 #include <trackbase/ActsGeometry.h>
0013
0014
0015 //TrkrCluster includes
0016 #include <trackbase/TrkrCluster.h>                      // for TrkrCluster
0017 #include <trackbase/TrkrDefs.h>                         // for getLayer, clu...
0018 #include <trackbase/TrkrClusterContainer.h>
0019 #include <trackbase/LaserClusterContainer.h>
0020
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/SubsysReco.h>                          // for SubsysReco
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHIODataNode.h>
0026 #include <phool/PHNode.h>                                // for PHNode
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/PHObject.h>                              // for PHObject
0029 #include <phool/PHRandomSeed.h>
0030 #include <phool/PHTimer.h>                               // for PHTimer
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>                                 // for PHWHERE
0033
0034 //ROOT includes for debugging
0035 #include <TFile.h>
0036 #include <TMatrixDSymfwd.h>                              // for TMatrixDSym
0037 #include <TMatrixTSym.h>                                 // for TMatrixTSym
0038 #include <TMatrixTUtils.h>                               // for TMatrixTRow
0039 #include <TNtuple.h>
0040 #include <TVector3.h>                                    // for TVector3
0041 #include <TVectorDfwd.h>                                 // for TVectorD
0042 #include <TVectorT.h>                                    // for TVectorT
0043
0044 // gsl
0045 #include <gsl/gsl_rng.h>
0046
0047 #include <Eigen/Core>
0048 #include <Eigen/Dense>
0049 #include <Eigen/LU>
0050
0051 //BOOST for combi seeding
0052 #include <boost/geometry/geometries/box.hpp>
0053 #include <boost/geometry/geometries/point.hpp>
0054 #include <boost/geometry/index/rtree.hpp>
0055
0056 // standard includes
0057 #include <algorithm>
0058 #include <assert.h>                                      // for assert
0059 #include <cfloat>
0060 #include <climits>                                      // for UINT_MAX
0061 #include <cmath>
0062 #include <cstdlib>                                      // for NULL, exit
0063 #include <fstream>
0064 #include <iostream>
0065 #include <iterator>                                      // for back_insert_...
0066 #include <memory>
0067 #include <tuple>
0068
0069 class TGeoManager;
0070
0071 using namespace Eigen;
0072
0073 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0074 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0075 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0076
0077
0078 using namespace std;
0079 namespace bg = boost::geometry;
0080 namespace bgi = boost::geometry::index;
0081 //end
0082
0083
0084
0085 #ifndef PHCOSMICSFILTER_H
0086 #define PHCOSMICSFILTER_H
0087
0088 #include "PHTrackSeeding.h"
0089
0090
0091 #if !defined(__CINT__) || defined(__CLING__)
0092 //BOOST for combi seeding
0093 #include <boost/geometry.hpp>
0094 #include <boost/geometry/geometries/box.hpp>
0095 #include <boost/geometry/geometries/point.hpp>
0096
0097 #include <boost/geometry/index/rtree.hpp>
0098 #endif
0099
0100
0101
0102
0103 // standard includes
0104 #include <cfloat>
0105 #include <iostream>                    // for operator<<, endl, basic_ostream
0106 #include <map>
0107 #include <string>                      // for string
0108 #include <vector>
0109
0110 // forward declarations
0111 class PHCompositeNode;
0112 class PHG4CellContainer;
0113 class PHG4CylinderGeomContainer;
0114 class PHG4HitContainer;
0115 class PHTimer;
0116 class sPHENIXSeedFinder;
0117 class SvtxClusterMap;
0118 class SvtxCluster;
0119 class SvtxTrackMap;
0120 class SvtxTrack;
0121 class SvtxVertexMap;
0122 class SvtxVertex;
0123 class TNtuple;
0124 class TFile;
0125 class TRKR_CLUSTER;
0126 class SvtxHitMap;
0127 class TrackSeedContainer;
0128
0129
0130 namespace bg = boost::geometry;
0131 namespace bgi = boost::geometry::index;
0132 typedef bg::model::point<float, 3, bg::cs::cartesian> point;
0133 typedef bg::model::box<point> box;
0134 typedef std::pair<point, int> pointKey;
0135 typedef std::pair<point, float> stubKey;
0136 typedef uint64_t cluskey;
0137
0138
0139
0140
0141 class PHStreakFinder : public SubsysReco
0142 {
0143  public:
0144   PHStreakFinder(const std::string &name = "PHStreakFinder",const std::string &filename = "streakntuple.root");
0145
0146   double chisq(const double *xx);
0147   //vector<TrkrCluster*> clusterpoints;
0148   virtual ~PHStreakFinder()
0149   {
0150   }
0151   void set_write_debug_ntuple(bool b){_write_ntp = b;}
0152   void set_create_tracks(bool b){_create_tracks = b;}
0153   void set_max_distance_to_origin(float val){ _max_dist_to_origin = val;}
0154   void set_min_nclusters(int n){ _min_nclusters = n;}
0155
0156  protected:
0157   int Setup(PHCompositeNode *topNode);
0158   int GetNodes(PHCompositeNode* topNode);
0159   int Init(PHCompositeNode *topNode) override;
0160   int InitRun(PHCompositeNode *topNode) override;
0161   int process_event(PHCompositeNode *topNode) override;
0162   int End(PHCompositeNode *topNode) override;
0163
0164
0165  private:
0166   /// fetch node pointers
0167   // int GetNodes(PHCompositeNode *topNode);
0168
0169    //static vector<tuple<double, double, double>> clusterpoints;
0170    /*static*/ //vector<TrkrCluster*> clusterpoints;
0171
0172   // node pointers
0173   LaserClusterContainer *_cluster_map = nullptr;
0174   //nodes to get norm vector
0175
0176   double phiadd(double phi1, double phi2);
0177   double phidiff(double phi1, double phi2);
0178   double pointKeyToTuple(pointKey *pK);
0179   double costfunction(const double *xx);
0180   //double chisq(const double *xx);
0181   void get_stub(const bgi::rtree<pointKey, bgi::quadratic<16>> &rtree, float pointx, float pointy, float pointz, int &count, double &xzslope, double &xzintercept, double &yzslope, double &yzintercept);
0182   ActsGeometry *tGeometry{nullptr};
0183 #ifndef __CINT__
0184  private:
0185   int createNodes(PHCompositeNode *topNode);
0186   std::string m_trackMapName = "TpcTrackSeedContainer";
0187   TrackSeedContainer *m_seedContainer = nullptr;
0188
0189   //int _nlayers_all;
0190   //unsigned int _nlayers_seeding;
0191   //std::vector<int> _seeding_layer;
0192   std::string _filename;
0193   unsigned int _nevent = 0;
0194   bool _write_ntp = false;
0195   bool _create_tracks = false;
0196   float _max_dist_to_origin = 0;
0197   unsigned int _min_nclusters = 50;
0198   TNtuple *_ntp_ev = nullptr;
0199   TNtuple *_ntp_clu = nullptr;
0200   TNtuple *_ntp_stub = nullptr;
0201   TNtuple *_ntp_trk = nullptr;
0202   TNtuple *_ntp_clutrk = nullptr;
0203   TFile *_tfile = nullptr;
0204   //std::vector<float> _radii_all;
0205
0206
0207 #endif  // __CINT__
0208 };
0209
0210 #endif