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 PHField;
0070 class TGeoManager;
0071
0072 using namespace Eigen;
0073
0074 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0075 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0076 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0077
0078
0079 using namespace std;
0080 namespace bg = boost::geometry;
0081 namespace bgi = boost::geometry::index;
0082 //end
0083
0084
0085
0086 #ifndef PHCOSMICSFILTER_H
0087 #define PHCOSMICSFILTER_H
0088
0089 #include "PHTrackSeeding.h"
0090
0091
0092 #if !defined(__CINT__) || defined(__CLING__)
0093 //BOOST for combi seeding
0094 #include <boost/geometry.hpp>
0095 #include <boost/geometry/geometries/box.hpp>
0096 #include <boost/geometry/geometries/point.hpp>
0097
0098 #include <boost/geometry/index/rtree.hpp>
0099 #endif
0100
0101
0102
0103
0104 // standard includes
0105 #include <cfloat>
0106 #include <iostream>                    // for operator<<, endl, basic_ostream
0107 #include <map>
0108 #include <string>                      // for string
0109 #include <vector>
0110
0111 // forward declarations
0112 class PHCompositeNode;
0113 class PHG4CellContainer;
0114 class PHG4CylinderGeomContainer;
0115 class PHG4HitContainer;
0116 class PHTimer;
0117 class sPHENIXSeedFinder;
0118 class SvtxClusterMap;
0119 class SvtxCluster;
0120 class SvtxTrackMap;
0121 class SvtxTrack;
0122 class SvtxVertexMap;
0123 class SvtxVertex;
0124 class TNtuple;
0125 class TFile;
0126 class TRKR_CLUSTER;
0127 class SvtxHitMap;
0128 class TrackSeedContainer;
0129
0130
0131 namespace bg = boost::geometry;
0132 namespace bgi = boost::geometry::index;
0133 typedef bg::model::point<float, 3, bg::cs::cartesian> point;
0134 typedef bg::model::box<point> box;
0135 typedef std::pair<point, int> pointKey;
0136 typedef std::pair<point, float> stubKey;
0137 typedef uint64_t cluskey;
0138
0139
0140
0141
0142 class PHStreakFinder : public SubsysReco
0143 {
0144  public:
0145   PHStreakFinder(const std::string &name = "PHStreakFinder",const std::string &filename = "streakntuple.root");
0146
0147   double chisq(const double *xx);
0148   //vector<TrkrCluster*> clusterpoints;
0149   virtual ~PHStreakFinder()
0150   {
0151   }
0152   void set_write_debug_ntuple(bool b){_write_ntp = b;}
0153   void set_create_tracks(bool b){_create_tracks = b;}
0154   void set_max_distance_to_origin(float val){ _max_dist_to_origin = val;}
0155   void set_min_nclusters(int n){ _min_nclusters = n;}
0156
0157  protected:
0158   int Setup(PHCompositeNode *topNode);
0159   int GetNodes(PHCompositeNode* topNode);
0160   int Init(PHCompositeNode *topNode) override;
0161   int InitRun(PHCompositeNode *topNode) override;
0162   int process_event(PHCompositeNode *topNode) override;
0163   int End(PHCompositeNode *topNode) override;
0164
0165
0166  private:
0167   /// fetch node pointers
0168   // int GetNodes(PHCompositeNode *topNode);
0169
0170    //static vector<tuple<double, double, double>> clusterpoints;
0171    /*static*/ //vector<TrkrCluster*> clusterpoints;
0172
0173   // node pointers
0174   LaserClusterContainer *_cluster_map = nullptr;
0175   //nodes to get norm vector
0176
0177   double phiadd(double phi1, double phi2);
0178   double phidiff(double phi1, double phi2);
0179   double pointKeyToTuple(pointKey *pK);
0180   double costfunction(const double *xx);
0181   //double chisq(const double *xx);
0182   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);
0183   ActsGeometry *tGeometry{nullptr};
0184 #ifndef __CINT__
0185  private:
0186   int createNodes(PHCompositeNode *topNode);
0187   std::string m_trackMapName = "TpcTrackSeedContainer";
0188   TrackSeedContainer *m_seedContainer = nullptr;
0189
0190   //int _nlayers_all;
0191   //unsigned int _nlayers_seeding;
0192   //std::vector<int> _seeding_layer;
0193   std::string _filename;
0194   unsigned int _nevent = 0;
0195   bool _write_ntp = false;
0196   bool _create_tracks = false;
0197   float _max_dist_to_origin = 0;
0198   unsigned int _min_nclusters = 50;
0199   TNtuple *_ntp_ev = nullptr;
0200   TNtuple *_ntp_clu = nullptr;
0201   TNtuple *_ntp_stub = nullptr;
0202   TNtuple *_ntp_trk = nullptr;
0203   TNtuple *_ntp_clutrk = nullptr;
0204   TFile *_tfile = nullptr;
0205   //std::vector<float> _radii_all;
0206
0207
0208 #endif  // __CINT__
0209 };
0210
0211 #endif