packages/trackreco/PHLineLaserReco.h

0001 #ifndef PHLINELASERRECO_H
0002 #define PHLINELASERRECO_H
0003
0004 /*!
0005  *  \file RPHLinelaserreco.h
0006  *  \RTree based hough tracking for cosmics
0007  *  \author Christof Roland
0008  */
0009
0010
0011 //begin
0012
0013 #include "PHTrackSeeding.h"
0014
0015 #include <fun4all/SubsysReco.h>
0016 #include <trackbase/TrkrDefs.h>  // for cluskey
0017 #include <trackbase/ActsGeometry.h>
0018
0019
0020 //TrkrCluster includes
0021 #include <trackbase/TrkrCluster.h>                      // for TrkrCluster
0022 #include <trackbase/TrkrDefs.h>                         // for getLayer, clu...
0023 #include <trackbase/TrkrClusterContainer.h>
0024 #include <trackbase/LaserClusterContainer.h>
0025
0026 #include <fun4all/Fun4AllReturnCodes.h>
0027 #include <fun4all/SubsysReco.h>                          // for SubsysReco
0028
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/PHIODataNode.h>
0031 #include <phool/PHNode.h>                                // for PHNode
0032 #include <phool/PHNodeIterator.h>
0033 #include <phool/PHObject.h>                              // for PHObject
0034 #include <phool/PHRandomSeed.h>
0035 #include <phool/PHTimer.h>                               // for PHTimer
0036 #include <phool/getClass.h>
0037 #include <phool/phool.h>                                 // for PHWHERE
0038
0039 //ROOT includes for debugging
0040 #include <TFile.h>
0041 #include <TMatrixDSymfwd.h>                              // for TMatrixDSym
0042 #include <TMatrixTSym.h>                                 // for TMatrixTSym
0043 #include <TMatrixTUtils.h>                               // for TMatrixTRow
0044 #include <TNtuple.h>
0045 #include <TVector3.h>                                    // for TVector3
0046 #include <TVectorDfwd.h>                                 // for TVectorD
0047 #include <TVectorT.h>                                    // for TVectorT
0048
0049 // gsl
0050 #include <gsl/gsl_rng.h>
0051
0052 #include <Eigen/Core>
0053 #include <Eigen/Dense>
0054 #include <Eigen/LU>
0055
0056 //BOOST for combi seeding
0057 #include <boost/geometry.hpp>
0058 #include <boost/geometry/geometries/box.hpp>
0059 #include <boost/geometry/geometries/point.hpp>
0060 #include <boost/geometry/index/rtree.hpp>
0061
0062
0063 // standard includes
0064 #include <algorithm>
0065 #include <cassert>                                      // for assert
0066 #include <cfloat>
0067 #include <climits>                                      // for UINT_MAX
0068 #include <cmath>
0069 #include <cstdlib>                                      // for NULL, exit
0070 #include <fstream>
0071 #include <iostream>
0072 #include <iterator>                                      // for back_insert_...
0073 #include <map>
0074 #include <memory>
0075 #include <string>                      // for string
0076 #include <tuple>
0077 #include <vector>
0078
0079
0080 class TGeoManager;
0081
0082
0083 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0084 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0085 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0086
0087
0088 using namespace Eigen;
0089 using namespace std;
0090 namespace bg = boost::geometry;
0091 namespace bgi = boost::geometry::index;
0092 //end
0093
0094 // forward declarations
0095 class PHCompositeNode;
0096 class PHG4CellContainer;
0097 class PHG4CylinderGeomContainer;
0098 class PHG4HitContainer;
0099 class PHTimer;
0100 class sPHENIXSeedFinder;
0101 class SvtxClusterMap;
0102 class SvtxCluster;
0103 class SvtxTrackMap;
0104 class SvtxTrack;
0105 class SvtxVertexMap;
0106 class SvtxVertex;
0107 class TNtuple;
0108 class TFile;
0109 class TRKR_CLUSTER;
0110 class SvtxHitMap;
0111 class TrackSeedContainer;
0112
0113
0114 namespace bg = boost::geometry;
0115 namespace bgi = boost::geometry::index;
0116 typedef bg::model::point<float, 3, bg::cs::cartesian> point;
0117 typedef bg::model::box<point> box;
0118 typedef std::pair<point, TrkrDefs::cluskey> pointInd;
0119
0120 typedef uint64_t cluskey;
0121
0122
0123
0124
0125 class PHLineLaserReco : public SubsysReco
0126 {
0127  public:
0128   PHLineLaserReco(const std::string &name = "PHLineLaserReco",
0129           const std::string &filename = "lasertuple.root"
0130           );
0131
0132   double chisq(const double *xx);
0133   //vector<TrkrCluster*> clusterpoints;
0134   virtual ~PHLineLaserReco()
0135   {
0136   }
0137   void set_write_debug_ntuple(bool b){_write_ntp = b;}
0138   void set_create_tracks(bool b){_create_tracks = b;}
0139   void set_max_distance_to_origin(float val){ _max_dist_to_origin = val;}
0140   void set_min_nclusters(int n){ _min_nclusters = n;}
0141   void segment(const int seg) { m_nseg = seg; }
0142   void runnumber(const int run) { m_nrun = run; }
0143   void job(const int job) { m_job = job; }
0144
0145  protected:
0146   int Setup(PHCompositeNode *topNode);
0147   int GetNodes(PHCompositeNode* topNode);
0148   int Init(PHCompositeNode *topNode) override;
0149   int InitRun(PHCompositeNode *topNode) override;
0150   int process_event(PHCompositeNode *topNode) override;
0151   int End(PHCompositeNode *topNode) override;
0152
0153
0154  private:
0155   int m_segment = 0;
0156   int m_runnumber = 0;
0157   int m_job = 0;
0158
0159   /// fetch node pointers
0160   // int GetNodes(PHCompositeNode *topNode);
0161
0162    //static vector<tuple<double, double, double>> clusterpoints;
0163    /*static*/ //vector<TrkrCluster*> clusterpoints;
0164
0165   // node pointers
0166   LaserClusterContainer *_cluster_map = nullptr;
0167   //nodes to get norm vector
0168
0169   double phiadd(double phi1, double phi2);
0170   double phidiff(double phi1, double phi2);
0171   double pointKeyToTuple(pointInd *pK);
0172   double costfunction(const double *xx);
0173   //double chisq(const double *xx);
0174   // void remove_stub(float msl, float mint, float rmsl2, float rmint2, std::multimap<unsigned short, stub> &all_stub_map);
0175   void get_stub(const bgi::rtree<pointInd, bgi::quadratic<16>> &rtree, float pointx, float pointy, int &count, double &slope, double &intercept);
0176   ActsGeometry *tGeometry{nullptr};
0177 #ifndef __CINT__
0178  private:
0179   int createNodes(PHCompositeNode *topNode);
0180   std::string m_trackMapName = "TpcTrackSeedContainer";
0181   TrackSeedContainer *m_seedContainer = nullptr;
0182
0183   //int _nlayers_all;
0184   //unsigned int _nlayers_seeding;
0185   //std::vector<int> _seeding_layer;
0186
0187   unsigned int _nevent = 0;
0188   bool _write_ntp = false;
0189   bool _create_tracks = false;
0190   float _max_dist_to_origin = 0;
0191   unsigned int _min_nclusters = 50;
0192   TNtuple *_ntp_cos = nullptr;
0193   TNtuple *_ntp_stub = nullptr;
0194   TNtuple *_ntp_max = nullptr;
0195   TNtuple *_ntp_trk = nullptr;
0196   TNtuple *_ntp_trk_hit = nullptr;
0197   TNtuple *_ntp_trk_clus = nullptr;
0198   TTree *m_hittree = nullptr;
0199   TTree *m_clustree = nullptr;
0200   TFile *_tfile = nullptr;
0201   std::string _filename;
0202   //std::vector<float> _radii_all;
0203   //! hit tree info
0204
0205   int m_nevent = std::numeric_limits<int>::quiet_NaN();
0206   int m_nrun = std::numeric_limits<int>::quiet_NaN();
0207   int m_nseg = std::numeric_limits<int>::quiet_NaN();
0208   int m_njob = std::numeric_limits<int>::quiet_NaN();
0209
0210   float m_hitx = std::numeric_limits<float>::quiet_NaN();
0211   float m_hity = std::numeric_limits<float>::quiet_NaN();
0212   float m_hitz = std::numeric_limits<float>::quiet_NaN();
0213   int m_hitadc = std::numeric_limits<int>::quiet_NaN();
0214   int m_hitlayer = std::numeric_limits<int>::quiet_NaN();
0215   int m_hitpad = std::numeric_limits<int>::quiet_NaN();
0216   int m_hittbin = std::numeric_limits<int>::quiet_NaN();
0217   int m_hitside = std::numeric_limits<int>::quiet_NaN();
0218   float m_slopexy = std::numeric_limits<float>::quiet_NaN();
0219   float m_interxy = std::numeric_limits<float>::quiet_NaN();
0220   float m_slopexz = std::numeric_limits<float>::quiet_NaN();
0221   float m_interxz = std::numeric_limits<float>::quiet_NaN();
0222   float m_slopeyz = std::numeric_limits<float>::quiet_NaN();
0223   float m_interyz = std::numeric_limits<float>::quiet_NaN();
0224   int m_nclus = std::numeric_limits<int>::quiet_NaN();
0225   float m_zfirst = std::numeric_limits<float>::quiet_NaN();
0226   float m_zlast = std::numeric_limits<float>::quiet_NaN();
0227
0228   float m_clux = std::numeric_limits<float>::quiet_NaN();
0229   float m_cluy = std::numeric_limits<float>::quiet_NaN();
0230   float m_cluz = std::numeric_limits<float>::quiet_NaN();
0231   int m_cluside = std::numeric_limits<int>::quiet_NaN();
0232   int m_cluadc = std::numeric_limits<int>::quiet_NaN();
0233   int m_clumaxadc = std::numeric_limits<int>::quiet_NaN();
0234   int m_size = std::numeric_limits<int>::quiet_NaN();
0235   int m_sizel = std::numeric_limits<int>::quiet_NaN();
0236   int m_sizephi = std::numeric_limits<int>::quiet_NaN();
0237   int m_sizet = std::numeric_limits<int>::quiet_NaN();
0238
0239 #endif  // __CINT__
0240 };
0241
0242 #endif