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 PHField;
0081 class TGeoManager;
0082
0083
0084 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0085 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0086 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0087
0088
0089 using namespace Eigen;
0090 using namespace std;
0091 namespace bg = boost::geometry;
0092 namespace bgi = boost::geometry::index;
0093 //end
0094
0095 // forward declarations
0096 class PHCompositeNode;
0097 class PHG4CellContainer;
0098 class PHG4CylinderGeomContainer;
0099 class PHG4HitContainer;
0100 class PHTimer;
0101 class sPHENIXSeedFinder;
0102 class SvtxClusterMap;
0103 class SvtxCluster;
0104 class SvtxTrackMap;
0105 class SvtxTrack;
0106 class SvtxVertexMap;
0107 class SvtxVertex;
0108 class TNtuple;
0109 class TFile;
0110 class TRKR_CLUSTER;
0111 class SvtxHitMap;
0112 class TrackSeedContainer;
0113
0114
0115 namespace bg = boost::geometry;
0116 namespace bgi = boost::geometry::index;
0117 typedef bg::model::point<float, 3, bg::cs::cartesian> point;
0118 typedef bg::model::box<point> box;
0119 typedef std::pair<point, TrkrDefs::cluskey> pointInd;
0120
0121 typedef uint64_t cluskey;
0122
0123
0124
0125
0126 class PHLineLaserReco : public SubsysReco
0127 {
0128  public:
0129   PHLineLaserReco(const std::string &name = "PHLineLaserReco",
0130           const std::string &filename = "lasertuple.root"
0131           );
0132
0133   double chisq(const double *xx);
0134   //vector<TrkrCluster*> clusterpoints;
0135   virtual ~PHLineLaserReco()
0136   {
0137   }
0138   void set_write_debug_ntuple(bool b){_write_ntp = b;}
0139   void set_create_tracks(bool b){_create_tracks = b;}
0140   void set_max_distance_to_origin(float val){ _max_dist_to_origin = val;}
0141   void set_min_nclusters(int n){ _min_nclusters = n;}
0142   void segment(const int seg) { m_nseg = seg; }
0143   void runnumber(const int run) { m_nrun = run; }
0144   void job(const int job) { m_job = job; }
0145
0146  protected:
0147   int Setup(PHCompositeNode *topNode);
0148   int GetNodes(PHCompositeNode* topNode);
0149   int Init(PHCompositeNode *topNode) override;
0150   int InitRun(PHCompositeNode *topNode) override;
0151   int process_event(PHCompositeNode *topNode) override;
0152   int End(PHCompositeNode *topNode) override;
0153
0154
0155  private:
0156   int m_segment = 0;
0157   int m_runnumber = 0;
0158   int m_job = 0;
0159
0160   /// fetch node pointers
0161   // int GetNodes(PHCompositeNode *topNode);
0162
0163    //static vector<tuple<double, double, double>> clusterpoints;
0164    /*static*/ //vector<TrkrCluster*> clusterpoints;
0165
0166   // node pointers
0167   LaserClusterContainer *_cluster_map = nullptr;
0168   //nodes to get norm vector
0169
0170   double phiadd(double phi1, double phi2);
0171   double phidiff(double phi1, double phi2);
0172   double pointKeyToTuple(pointInd *pK);
0173   double costfunction(const double *xx);
0174   //double chisq(const double *xx);
0175   // void remove_stub(float msl, float mint, float rmsl2, float rmint2, std::multimap<unsigned short, stub> &all_stub_map);
0176   void get_stub(const bgi::rtree<pointInd, bgi::quadratic<16>> &rtree, float pointx, float pointy, int &count, double &slope, double &intercept);
0177   ActsGeometry *tGeometry{nullptr};
0178 #ifndef __CINT__
0179  private:
0180   int createNodes(PHCompositeNode *topNode);
0181   std::string m_trackMapName = "TpcTrackSeedContainer";
0182   TrackSeedContainer *m_seedContainer = nullptr;
0183
0184   //int _nlayers_all;
0185   //unsigned int _nlayers_seeding;
0186   //std::vector<int> _seeding_layer;
0187
0188   unsigned int _nevent = 0;
0189   bool _write_ntp = false;
0190   bool _create_tracks = false;
0191   float _max_dist_to_origin = 0;
0192   unsigned int _min_nclusters = 50;
0193   TNtuple *_ntp_cos = nullptr;
0194   TNtuple *_ntp_stub = nullptr;
0195   TNtuple *_ntp_max = nullptr;
0196   TNtuple *_ntp_trk = nullptr;
0197   TNtuple *_ntp_trk_hit = nullptr;
0198   TNtuple *_ntp_trk_clus = nullptr;
0199   TTree *m_hittree = nullptr;
0200   TTree *m_clustree = nullptr;
0201   TFile *_tfile = nullptr;
0202   std::string _filename;
0203   //std::vector<float> _radii_all;
0204   //! hit tree info
0205
0206   int m_nevent = std::numeric_limits<int>::quiet_NaN();
0207   int m_nrun = std::numeric_limits<int>::quiet_NaN();
0208   int m_nseg = std::numeric_limits<int>::quiet_NaN();
0209   int m_njob = std::numeric_limits<int>::quiet_NaN();
0210
0211   float m_hitx = std::numeric_limits<float>::quiet_NaN();
0212   float m_hity = std::numeric_limits<float>::quiet_NaN();
0213   float m_hitz = std::numeric_limits<float>::quiet_NaN();
0214   int m_hitadc = std::numeric_limits<int>::quiet_NaN();
0215   int m_hitlayer = std::numeric_limits<int>::quiet_NaN();
0216   int m_hitpad = std::numeric_limits<int>::quiet_NaN();
0217   int m_hittbin = std::numeric_limits<int>::quiet_NaN();
0218   int m_hitside = std::numeric_limits<int>::quiet_NaN();
0219   float m_slopexy = std::numeric_limits<float>::quiet_NaN();
0220   float m_interxy = std::numeric_limits<float>::quiet_NaN();
0221   float m_slopexz = std::numeric_limits<float>::quiet_NaN();
0222   float m_interxz = std::numeric_limits<float>::quiet_NaN();
0223   float m_slopeyz = std::numeric_limits<float>::quiet_NaN();
0224   float m_interyz = std::numeric_limits<float>::quiet_NaN();
0225   int m_nclus = std::numeric_limits<int>::quiet_NaN();
0226   float m_zfirst = std::numeric_limits<float>::quiet_NaN();
0227   float m_zlast = std::numeric_limits<float>::quiet_NaN();
0228
0229   float m_clux = std::numeric_limits<float>::quiet_NaN();
0230   float m_cluy = std::numeric_limits<float>::quiet_NaN();
0231   float m_cluz = std::numeric_limits<float>::quiet_NaN();
0232   int m_cluside = std::numeric_limits<int>::quiet_NaN();
0233   int m_cluadc = std::numeric_limits<int>::quiet_NaN();
0234   int m_clumaxadc = std::numeric_limits<int>::quiet_NaN();
0235   int m_size = std::numeric_limits<int>::quiet_NaN();
0236   int m_sizel = std::numeric_limits<int>::quiet_NaN();
0237   int m_sizephi = std::numeric_limits<int>::quiet_NaN();
0238   int m_sizet = std::numeric_limits<int>::quiet_NaN();
0239
0240 #endif  // __CINT__
0241 };
0242
0243 #endif