Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:30 

0001 /*!
0002  *  \file PHLineLaserReco.cc
0003  
0004  *  \author Christof Roland 
0005  */
0006 
0007 
0008 //begin
0009 
0010 #include "PHLineLaserReco.h"
0011 
0012 #include "AssocInfoContainer.h"                         // for AssocInfoCont...
0013 
0014 // sPHENIX Geant4 includes
0015 #include <g4detectors/PHG4CylinderCellGeom.h>
0016 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0017 #include <g4detectors/PHG4CylinderGeom.h>
0018 #include <g4detectors/PHG4CylinderGeomContainer.h>
0019 
0020 #include <trackbase/TrackFitUtils.h>
0021 #include <trackbase/TrkrCluster.h>
0022 #include <trackbase/TrkrClusterContainer.h>
0023 #include <trackbase/TrkrClusterHitAssoc.h>
0024 #include <trackbase/TrkrClusterIterationMapv1.h>
0025 #include <trackbase/TrkrDefs.h>  // for getLayer, clu...
0026 #include <trackbase/TpcDefs.h>  // for getLayer, clu...
0027 #include <trackbase/LaserCluster.h>
0028 #include <trackbase/LaserClusterContainer.h>
0029 #include <trackbase/LaserClusterContainerv1.h>
0030 #include <trackbase/LaserClusterv1.h>
0031 
0032 #include <trackbase_historic/TrackSeed.h>
0033 #include <trackbase_historic/TrackSeedContainer.h>
0034 #include <trackbase_historic/TrackSeedContainer_v1.h>
0035 #include <trackbase_historic/TrackSeed_v2.h>
0036 
0037 // sPHENIX includes
0038 #include <fun4all/Fun4AllReturnCodes.h>
0039 #include <fun4all/SubsysReco.h>
0040 
0041 #include <phool/PHTimer.h>                              // for PHTimer
0042 #include <phool/getClass.h>
0043 
0044 //ROOT includes for debugging
0045 #include <TFile.h>
0046 #include <TH1.h>
0047 #include <TNtuple.h>
0048 #include <TAxis.h>
0049 #include <TGraph.h>
0050 
0051 #include <Eigen/Core>                  // for Matrix
0052 #include <Eigen/Dense>
0053 
0054 //BOOST for combi seeding
0055 #include <boost/geometry.hpp>
0056 #include <boost/geometry/geometries/box.hpp>
0057 #include <boost/geometry/geometries/point.hpp>
0058 #include <boost/geometry/index/rtree.hpp>
0059 
0060 // standard includes
0061 #include <algorithm>
0062 #include <cfloat>
0063 #include <climits>                                     // for UINT_MAX
0064 #include <cmath>
0065 #include <iostream>
0066 #include <list>
0067 #include <map>
0068 #include <memory>
0069 #include <memory>
0070 #include <set>                                          // for set
0071 #include <string>                                  // for string
0072 #include <tuple>
0073 #include <utility>
0074 #include <vector>
0075                                    // for pair, make_pair
0076 
0077 // forward declarations
0078 class BbcVertexMap;
0079 
0080 class PHCompositeNode;
0081 
0082 class PHG4CellContainer;
0083 class PHG4CylinderGeomContainer;
0084 class PHG4Particle;
0085 class PHG4TruthInfoContainer;
0086 class PHTimer;
0087 
0088 class sPHENIXSeedFinder;
0089 
0090 class SvtxClusterMap;
0091 class SvtxCluster;
0092 class SvtxTrackMap;
0093 class SvtxTrack;
0094 class SvtxTrackState;
0095 class SvtxVertexMap;
0096 class SvtxHitMap;
0097 
0098 class TNtuple;
0099 class TFile;
0100 
0101 
0102 using point = bg::model::point<float, 3, bg::cs::cartesian>;
0103 //typedef bg::model::point<float, 3, bg::cs::cartesian> point;
0104 using box = bg::model::box<point>;
0105 //typedef bg::model::box<point> box;
0106 using pointInd = std::pair<point, TrkrDefs::cluskey>;
0107 //typedef std::pair<point, TrkrDefs::cluskey> pointKey;
0108 
0109 using myrtree = bgi::rtree<pointInd, bgi::quadratic<16>>;
0110 
0111 
0112 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0113 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0114 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0115 
0116 
0117 using namespace std;
0118 
0119 vector<LaserCluster*> laserpoints;
0120 
0121 PHLineLaserReco::PHLineLaserReco(const std::string& name, const string& filename)
0122   : SubsysReco(name)
0123   , _filename(filename)
0124 {
0125   
0126 }
0127 
0128 int PHLineLaserReco::GetNodes(PHCompositeNode* topNode)
0129 {
0130   tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0131   if (!tGeometry)
0132   {
0133     std::cout << PHWHERE << "No acts reco geometry, bailing."
0134               << std::endl;
0135     return Fun4AllReturnCodes::ABORTEVENT;
0136   }
0137   _cluster_map = findNode::getClass<LaserClusterContainerv1>(topNode, "LASER_CLUSTER");
0138   if (!_cluster_map)
0139   {
0140     std::cout << PHWHERE << "No cluster container on node tree, bailing."
0141               << std::endl;
0142     return Fun4AllReturnCodes::ABORTEVENT;
0143   }
0144   return Fun4AllReturnCodes::EVENT_OK;
0145 }
0146 
0147 //____________________________________________________________________________..
0148 int PHLineLaserReco::Init(PHCompositeNode* topNode)
0149 {
0150   if(topNode==nullptr){
0151     std::cout << PHWHERE << "No topNode, bailing."
0152           << std::endl;
0153     return Fun4AllReturnCodes::ABORTEVENT;
0154   }
0155   _nevent = 0;
0156   if(_write_ntp){
0157   _tfile = new TFile(_filename.c_str(), "RECREATE");
0158   _ntp_cos = new TNtuple("ntp_cos", "cos event info","ev:x:y:z:adc:maxadc:size:lsize:phisize:tsize:side");
0159   _ntp_stub = new TNtuple("ntp_stub", "cos stub info","ev:int:sl:tan:x0:count");
0160   _ntp_max = new TNtuple("ntp_max", "cos stub info","ev:intmin:intmax:slmin:slmax");
0161   _ntp_trk = new TNtuple("ntp_trk", "laser track info","ev:int:sl:intfit:slfit:nclus");
0162   //  _ntp_trk_clus = new TNtuple("ntp_trk_clus", "laser clus","ev:x:y:z:adc:maxadc:size:lsize:phisize:tsize:intfit:slfit:nclus:zfirst:zlast");
0163   _ntp_trk_hit = new TNtuple("ntp_trk_hit", "laser hitinfo","ev:x:y:z:adc:lay:phibin:zbin:intfit:slfit:nclus:zfirst:zlast");
0164   }
0165 
0166   m_hittree = new TTree("hittree", "A tree with all hits");
0167   m_hittree->Branch("event", &m_nevent, "m_nevent/I");
0168   m_hittree->Branch("run", &m_nrun, "m_nrun/I");
0169   m_hittree->Branch("seg", &m_nseg, "m_nseg/I");
0170   m_hittree->Branch("job", &m_njob, "m_njob/I");
0171   m_hittree->Branch("x", &m_hitx, "m_hitx/F");
0172   m_hittree->Branch("y", &m_hity, "m_hity/F");
0173   m_hittree->Branch("z", &m_hitz, "m_hitz/F");
0174   m_hittree->Branch("layer", &m_hitlayer, "m_hitlayer/I");
0175   m_hittree->Branch("adc", &m_hitadc, "m_hitadc/I");
0176   m_hittree->Branch("pad", &m_hitpad, "m_hitpad/I");
0177   m_hittree->Branch("tbin", &m_hittbin, "m_hittbin/I");
0178   m_hittree->Branch("side", &m_hitside, "m_hitside/I");
0179   m_hittree->Branch("slopexy", &m_slopexy, "m_slopexy/F");
0180   m_hittree->Branch("interxy", &m_interxy, "m_interxy/F");
0181   m_hittree->Branch("slopexz", &m_slopexz, "m_slopexz/F");
0182   m_hittree->Branch("interxz", &m_interxz, "m_interxz/F");
0183   m_hittree->Branch("slopeyz", &m_slopeyz, "m_slopeyz/F");
0184   m_hittree->Branch("interyz", &m_interyz, "m_interyz/F");
0185   m_hittree->Branch("nclus", &m_nclus, "m_nclus/I");
0186   m_hittree->Branch("zlast", &m_zlast, "m_zlast/F");
0187   m_hittree->Branch("zfirst", &m_zfirst, "m_first/F");
0188 
0189   m_clustree = new TTree("clustree", "A tree with all hits");
0190   m_clustree->Branch("event", &m_nevent, "m_nevent/I");
0191   m_clustree->Branch("run", &m_nrun, "m_nrun/I");
0192   m_clustree->Branch("seg", &m_nseg, "m_nseg/I");
0193   m_clustree->Branch("job", &m_njob, "m_njob/I");
0194   m_clustree->Branch("x", &m_clux, "m_clux/F");
0195   m_clustree->Branch("y", &m_cluy, "m_cluy/F");
0196   m_clustree->Branch("z", &m_cluz, "m_cluz/F");
0197   m_clustree->Branch("adc", &m_cluadc, "m_cluadc/I");
0198   m_clustree->Branch("side", &m_cluside, "m_cluside/I");
0199   m_clustree->Branch("maxadc", &m_clumaxadc, "m_clumaxadc/I");
0200   m_clustree->Branch("size", &m_size, "m_size/I");
0201   m_clustree->Branch("lsize", &m_sizel, "m_sizel/I");
0202   m_clustree->Branch("phisize", &m_sizephi, "m_sizephi/I");
0203   m_clustree->Branch("tsize", &m_sizet, "m_sizet/I");
0204   m_clustree->Branch("slopexy", &m_slopexy, "m_slopexy/F");
0205   m_clustree->Branch("interxy", &m_interxy, "m_interxy/F");
0206   m_clustree->Branch("slopexz", &m_slopexz, "m_slopexz/F");
0207   m_clustree->Branch("interxz", &m_interxz, "m_interxz/F");
0208   m_clustree->Branch("slopeyz", &m_slopeyz, "m_slopeyz/F");
0209   m_clustree->Branch("interyz", &m_interyz, "m_interyz/F");
0210   m_clustree->Branch("nclus", &m_nclus, "m_nclus/I");
0211   m_clustree->Branch("zlast", &m_zlast, "m_zlast/F");
0212   m_clustree->Branch("zfirst", &m_zfirst, "m_first/F");
0213   
0214   return Fun4AllReturnCodes::EVENT_OK;
0215 }
0216 
0217 //____________________________________________________________________________..
0218 int PHLineLaserReco::InitRun(PHCompositeNode* topNode)
0219 {
0220   int ret = GetNodes(topNode);
0221   if (ret != Fun4AllReturnCodes::EVENT_OK)
0222   {
0223     return ret;
0224   }
0225   ret = createNodes(topNode);
0226 
0227 
0228   return ret;
0229 }
0230 int PHLineLaserReco::createNodes(PHCompositeNode* topNode)
0231 {
0232   PHNodeIterator iter(topNode);
0233 
0234   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0235 
0236   if (!dstNode)
0237   {
0238     std::cerr << "DST node is missing, quitting" << std::endl;
0239     throw std::runtime_error("Failed to find DST node in PHLineLaserReco::createNodes");
0240   }
0241 
0242   PHNodeIterator dstIter(dstNode);
0243   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0244 
0245   if (!svtxNode)
0246   {
0247     svtxNode = new PHCompositeNode("SVTX");
0248     dstNode->addNode(svtxNode);
0249   }
0250 
0251   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
0252   if (!m_seedContainer)
0253   {
0254     m_seedContainer = new TrackSeedContainer_v1;
0255     PHIODataNode<PHObject>* trackNode =
0256         new PHIODataNode<PHObject>(m_seedContainer, m_trackMapName, "PHObject");
0257     svtxNode->addNode(trackNode);
0258   }
0259   if (m_seedContainer){
0260     if(Verbosity()>0){ cout << "SEED CONTAINER CREATED" << endl;}
0261   }
0262   return Fun4AllReturnCodes::EVENT_OK;
0263 }
0264 
0265 double PHLineLaserReco::phiadd(double phi1, double phi2){
0266   double s = phi1+phi2;
0267   if(s>2*M_PI) {return s-2*M_PI;}
0268   else if(s<0) {return s+2*M_PI;}
0269   else {return s;}
0270 }
0271 
0272 
0273 double PHLineLaserReco::phidiff(double phi1, double phi2){
0274   double d = phi1-phi2;
0275   if(d>M_PI) {return d-2*M_PI;}
0276   else if(d<-M_PI) {return d+2*M_PI;}
0277   else {return d;}
0278 }
0279 
0280 void PHLineLaserReco::get_stub(const myrtree &search_rtree, float pointx, float pointy, int &count, double &slope, double &intercept){ //NOLINT
0281   float m1_dx = 6;
0282   float m1_dy = 6;
0283   vector<pointInd> boxclusters;
0284   search_rtree.query(bgi::intersects(box(point(pointx-m1_dx,pointy-m1_dy,-1),
0285                   point(pointx+m1_dx,pointy+m1_dy,1))),std::back_inserter(boxclusters));
0286   
0287   int nbox = boxclusters.size();
0288   
0289   int nhit = 0;
0290   double xsum = 0;
0291   double x2sum = 0;
0292   double ysum = 0;
0293   double xysum = 0;
0294   if(nbox>=2){
0295     for(vector<pointInd>::iterator pbox = boxclusters.begin();pbox!=boxclusters.end();++pbox){
0296       float boxx = pbox->first.get<0>();
0297       float boxy = pbox->first.get<1>();
0298       //      std::cout << "fprobe: " << count << " x:"  << boxx << " y: " << boxy << std::endl;
0299       nhit++;
0300       float r = boxx;
0301       float z = boxy;
0302       xsum = xsum + r;            // calculate sigma(xi)
0303       ysum = ysum + z;            // calculate sigma(yi)
0304       x2sum = x2sum + r*r;  // calculate sigma(x^2i)
0305       xysum = xysum + r * z;      // calculate sigma(xi*yi)
0306     }
0307   }
0308   
0309   const double denominator = ((x2sum * nhit) - (xsum*xsum));
0310   slope = (xysum * nhit - xsum * ysum) / denominator;   // calculate slope
0311   intercept = (x2sum * ysum - xsum * xysum) / denominator;  // calculate intercept
0312   count = nhit;
0313 }
0314 
0315 int PHLineLaserReco::process_event(PHCompositeNode* topNode)
0316 {
0317   if(Verbosity()>0){std::cout << "start processing event: " << std::endl;}
0318   
0319   if(topNode==nullptr){
0320     std::cout << PHWHERE << "No topNode, bailing."
0321           << std::endl;
0322     return Fun4AllReturnCodes::ABORTEVENT;
0323   }
0324     
0325   _nevent++;
0326   std::vector<TrackSeed_v2> clean_chains;
0327   //Fill rtree
0328   bgi::rtree<pointInd, bgi::quadratic<16> > rtree;
0329 
0330   float slmin =  99999999999.9;
0331   float slmax = -99999999999.9;
0332   float intmin = 99999999999.9;
0333   float intmax =-99999999999.9;
0334 //  int num = 0;
0335   //  for(const auto& hitsetkey:_cluster_map->getHitSetKeys(TrkrDefs::TrkrId::tpcId)){
0336   if(_cluster_map->size()<_min_nclusters){
0337     if(Verbosity()>0){std::cout << " not enough clusters in event: " << _cluster_map->size() << " < " << _min_nclusters << endl;}
0338     return Fun4AllReturnCodes::ABORTEVENT;
0339   }
0340 
0341   int numberofseeds = 0;
0342 
0343   auto range = _cluster_map->getClusters();
0344   if(Verbosity()>0){std::cout << "_nevent " << _nevent << "n clusters" << std::distance(range.first,range.second) << std::endl;}
0345   for(unsigned int process_side=0;process_side<2;process_side++){
0346   for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
0347     LaserCluster *cluster = clusIter->second;
0348     TrkrDefs::cluskey lckey = clusIter->first;
0349     unsigned int tpcside = TpcDefs::getSide(lckey);
0350     if(tpcside!=process_side){
0351       continue;
0352     }
0353     
0354     const Acts::Vector3 globalpos_d =  {cluster->getX() , cluster->getY(), cluster->getZ()};
0355     const Acts::Vector3 globalpos = { globalpos_d.x(), globalpos_d.y(), globalpos_d.z()};
0356     
0357     //    int side = -1;
0358     unsigned int adcSum = 0.0;
0359     unsigned int maxAdc = 0.0;
0360     int iphimin = 6666, iphimax = -1;
0361     int ilaymin = 6666, ilaymax = -1;
0362     float itmin = 66666666.6, itmax = -6666666666.6;
0363     unsigned int nHits = cluster->getNhits();
0364     
0365     for(unsigned int i  = 0; i< nHits;i++){
0366       //      float coords[3] = {cluster->getHitX(i), cluster->getHitY(i), cluster->getHitZ(i)};
0367       int iphi = cluster->getHitIPhi(i);
0368       int it = cluster->getHitIT(i);
0369       int ilay = cluster->getHitLayer(i);
0370       unsigned int adc =  cluster->getHitAdc(i);
0371       if(iphi<iphimin){iphimin = iphi;}
0372       if(iphi>iphimax){iphimax = iphi;}
0373       if(ilay<ilaymin){ilaymin = ilay;}
0374       if(ilay>ilaymax){ilaymax = ilay;}
0375       if(it<itmin){itmin = it;}
0376       if(it>itmax){itmax = it;}
0377       adcSum += adc;
0378       if(adc>maxAdc){
0379     maxAdc = adc;
0380       }
0381     }
0382     int phisize =  iphimax - iphimin + 1;
0383     int lsize =  ilaymax - ilaymin + 1;
0384     int tsize = itmax - itmin +1;
0385     
0386     if(_write_ntp)
0387       //    ev:x:y:z:adc:maxadc:size:lsize:phisize:tsize
0388       _ntp_cos->Fill(_nevent,globalpos_d.x(), globalpos_d.y(),globalpos_d.z(),adcSum,maxAdc,nHits,lsize,phisize,tsize,tpcside);
0389     if(lsize>2||maxAdc<50||tsize<=1)
0390       continue;
0391     std::vector<pointInd> testduplicate;
0392     rtree.query(bgi::intersects(box(point(globalpos_d.x()-0.001,globalpos_d.y()-0.001,-1),
0393                     point(globalpos_d.x()+0.001,globalpos_d.y()+0.001,1))),std::back_inserter(testduplicate));
0394     if (!testduplicate.empty()){
0395       continue;
0396     }
0397     rtree.insert(std::make_pair(point(globalpos.x() , globalpos.y(), 0.0), lckey)); 
0398   }
0399   
0400   //Get all clusters from rtree, fit stubs around clusters, fill stub tree
0401   vector<pointInd> allclusters;
0402   
0403   rtree.query(bgi::intersects(box(point(-80,-80,-1),point(80,80,1))),std::back_inserter(allclusters));
0404   if(Verbosity()>0) cout << "number clus is " << allclusters.size() << endl;
0405   
0406   bgi::rtree<pointInd, bgi::quadratic<16> > rtree_stub;
0407   
0408   for(vector<pointInd>::iterator cluster = allclusters.begin();cluster!=allclusters.end();++cluster){
0409     float pointx = cluster->first.get<0>();
0410     float pointy = cluster->first.get<1>();
0411     int fcount = 0;
0412     double fslope = 0;
0413     double fintercept = 0;
0414     //calc slope and intersect from 5 cluster stubs
0415     get_stub(rtree, pointx, pointy, fcount, fslope, fintercept);
0416     if(finite(fslope)&&fcount>5){
0417       rtree_stub.insert(std::make_pair(point(fintercept , fslope, 0.0), 0));
0418       float tana = atan(fslope);
0419       if(_write_ntp)
0420     _ntp_stub->Fill(_nevent,fintercept,fslope,tana,fintercept,fcount);
0421       if(fslope>slmax)slmax = fslope;
0422       if(fslope<slmin)slmin = fslope;
0423       if(fintercept>intmax)intmax = fintercept;
0424       if(fintercept<intmin)intmin = fintercept;
0425     }
0426   }
0427   //  int out = 0;
0428   // find clusters of slope intectept pairs
0429   map <int,pair<float,float>> outtrkmap;
0430   int nout = 0;
0431   
0432   while(rtree_stub.size()>3){
0433     std::vector<pointInd> allstubs;
0434     rtree_stub.query(bgi::intersects(box(point(intmin,slmin,-1),point(intmax,slmax,1))),std::back_inserter(allstubs));
0435     map <int,pair<float,float>> trkmap;
0436     
0437     float int_width = (intmax - intmin)/10;
0438     float sl_width = (slmax - slmin)/10;
0439     for(vector<pointInd>::iterator stub = allstubs.begin();stub!=allstubs.end();++stub){
0440       float pint = stub->first.get<0>();
0441       float psl = stub->first.get<1>();
0442       vector<pointInd> trkcand;
0443       
0444       rtree_stub.query(bgi::intersects(box(point(pint-int_width,psl-sl_width,-1),point(pint+int_width,psl+sl_width,1))),std::back_inserter(trkcand));
0445       
0446       int ntrk = trkcand.size();
0447       int count = 0;
0448       float intsum = 0;
0449       float slsum = 0;
0450       if(ntrk>=5){
0451     for(vector<pointInd>::iterator ptrk = trkcand.begin();ptrk!=trkcand.end();++ptrk){
0452       float trkint = ptrk->first.get<0>();
0453       float trksl = ptrk->first.get<1>();
0454       if(Verbosity()>0) cout<< "    stub " << ntrk
0455                 << " int: " << trkint 
0456                 << " sl: "  << trksl 
0457                 << endl;
0458       
0459       intsum = intsum + trkint;            // calculate sigma(xi)
0460       slsum = slsum + trksl;            // calculate sigma(yi)
0461       count++;
0462     }
0463     float mint = (intsum/count);
0464     float msl  = (slsum/count);
0465     trkmap[ntrk] = std::make_pair(mint,msl); 
0466     
0467     if(Verbosity()>0) cout<< " stub in box " << ntrk
0468                   << " int: " << pint 
0469                   << " sl: " << psl 
0470                   << " intwi: " << int_width 
0471                   << " slwi: " << sl_width
0472                   << " mint " << mint
0473                   << " msl " << msl
0474                   << endl;
0475       }
0476     }
0477     
0478     if( trkmap.size()>0){
0479       int size_before = rtree_stub.size();
0480       if(Verbosity()>0) cout << "mapend: " << trkmap.rbegin()->first << " | " << (trkmap.rbegin()->second).first << " | "  << trkmap.rbegin()->second.second << endl;
0481       // for (const auto& [key, value] : trkmap)
0482       //    std::cout <<" ev: " << _nevent << '[' << key << "] = " << value.first << " | " << value.second << std::endl;
0483       //remove stubs in box
0484       float rmint = trkmap.rbegin()->second.first;
0485       float rmsl = trkmap.rbegin()->second.second;
0486       vector<pointInd> rmcand;
0487       rtree_stub.query(bgi::intersects(box(point(rmint-int_width,rmsl-sl_width,-1),point(rmint+int_width,rmsl+sl_width,1))),std::back_inserter(rmcand));
0488       for(vector<pointInd>::iterator rmstub = rmcand.begin();rmstub!=rmcand.end();++rmstub){
0489     float rmpint = rmstub->first.get<0>();
0490     float rmpsl = rmstub->first.get<1>();
0491     if(Verbosity()>0) cout<< "    rm " <<  " int: " << rmpint 
0492                   << " sl: " << rmpsl 
0493                   << endl;
0494     
0495     rtree_stub.remove(*rmstub);
0496     /*  if (!rtree_stub.query(bgi::nearest(rmpoint, 1), &match))
0497         continue;
0498     */
0499     /*
0500       if (bg::distance(rmpoint, match.first) <= 0.5) {
0501       rtree_stub.remove(match);
0502       
0503     */
0504       }
0505       int size_after = rtree_stub.size();
0506       if(size_before == size_after) {break;}
0507       outtrkmap.insert(std::make_pair(nout,std::make_pair(rmint,rmsl)));
0508       nout++;
0509       if(Verbosity()>0){cout << " tree sixe after remove: " << rtree_stub.size() << endl;}
0510     }else{
0511       break;
0512     }
0513     trkmap.clear();
0514   }
0515   
0516   //  bool keep_event = false;
0517   for (const auto& [key, value] : outtrkmap){
0518     if(Verbosity()>0) std::cout <<" out ev: " << _nevent << '[' << key << "] = " << value.first << " | " << value.second << std::endl;
0519     //Loop over clusters and pick up clusters compatible with line parameters
0520     
0521     float tint = -1*value.first;
0522     float tsl = -1*value.second;
0523     
0524     //  float liney = tsl*ptx + tint;
0525     
0526     // double ry = pty;
0527     // double zx = ptx;
0528     
0529     // The shortest distance of a point from a circle is along the radial; line from the circle center to the point
0530     
0531     double a = -tsl;//-slope;
0532     double b = -1.0;
0533     double c = -tint;//-intercept;
0534     double r = 100;
0535     double x0 = -a*c/(a*a+b*b), y0 = -b*c/(a*a+b*b);
0536     double ax = NAN, ay = NAN, bx = NAN, by = NAN;
0537     if ((c*c) > (r*r*(a*a+b*b)+0.00000001)){
0538       if(Verbosity()>0) cout << "no points " << endl;
0539     }
0540     else if (abs (c*c - r*r*(a*a+b*b)) < 0.0000001) {
0541       if(Verbosity()>0){ puts ("1 point");
0542     cout << x0 << ' ' << y0 << endl;
0543       }
0544     }
0545     else {
0546       double d = r*r - c*c/(a*a+b*b);
0547       double mult = sqrt (d / (a*a+b*b));
0548       ax = x0 + b * mult;
0549       bx = x0 - b * mult;
0550       ay = y0 - a * mult;
0551       by = y0 + a * mult;
0552       if(Verbosity()>0){puts ("2 points");
0553     cout << ax << ' ' << ay << '\n' << bx << ' ' << by << endl;
0554       }
0555     }
0556     float dist = 4;
0557     float alpha = atan(tsl);
0558     float offax = -1*dist*sin(alpha);
0559     float offay = -1*dist*cos(alpha);
0560     float offbx = 1*dist*sin(alpha);
0561     float offby = 1*dist*cos(alpha);
0562     
0563     //float trkdist = 4;
0564     vector<pointInd> lineclusters;
0565     if(Verbosity()>0){ cout << ax << ' ' << ay << ' ' << bx << ' ' << by << endl;
0566       cout << "auto b0 = new TLine(" << ax << ',' << ay << ',' << bx  << ',' << by <<")" << endl;
0567       cout << "auto b1 = new TLine(" << ax +offax << ',' << ay +offay << ',' << bx -offbx << ',' << by - offby<<")" << endl;
0568       cout << "auto b2 = new TLine(" << ax -offax << ',' << ay -offay << ',' << bx +offbx << ',' << by + offby<<")" << endl;
0569       cout << "b0->Draw(\"same\")" << endl;
0570       cout << "b1->Draw(\"same\")" << endl;
0571       cout << "b2->Draw(\"same\")" << endl;
0572     }
0573     rtree.query(bgi::intersects(box(point(-80,-80,-1),point(80,80,1))),std::back_inserter(lineclusters));
0574     if(Verbosity()>0)  cout << "number line clus is " << lineclusters.size() << endl;
0575     
0576     // Check if track hits MVTX (dist to 0,0 < 2)
0577     
0578     //first fit
0579     int nhitff = 0;
0580     double xsumff = 0;
0581     double ysumff = 0;
0582     double zsumff = 0;
0583     double x2sumff = 0;
0584     double y2sumff = 0;
0585     double z2sumff = 0;
0586     double xysumff = 0;
0587     double xzsumff = 0;
0588     double yzsumff = 0;
0589     double zfirst = 66666666;
0590     double zlast = -66666666;
0591     vector<pointInd> trkclusters;
0592     unsigned int trk_nclus = 0;
0593     
0594     //    fitstraight lines first
0595     for(vector<pointInd>::iterator clustertrk = lineclusters.begin();clustertrk!=lineclusters.end();++clustertrk){
0596       
0597       float ptx = clustertrk->first.get<0>();
0598       float pty = clustertrk->first.get<1>();
0599       TrkrDefs::cluskey tcind = clustertrk->second;
0600       LaserCluster *las_clus = _cluster_map->findCluster(tcind);      
0601       float ptz = las_clus->getZ();
0602       
0603       float res =  std::abs(a*ptx+b*pty+c)/sqrt(a*a+b*b);
0604       //float res = std::abs(tsl*ptx+(-1*pty)+tint)/sqrt((tsl*tsl+1));
0605       // if(Verbosity()>0)
0606       if(Verbosity()>0){cout << " x: " << ptx << " y: " << pty << " res " << res << endl;}
0607       if(res<4){
0608     float boxx = ptx;
0609     float boxy = pty;
0610     float boxz = ptz;
0611     if(Verbosity()>0){std::cout << "fprobe x:"  << boxx << " y: " << boxy <<  " z: " << boxz << std::endl;}
0612     nhitff++;
0613     float xff = boxx;
0614     float yff = boxy;
0615     float zff = boxz;
0616     
0617     xsumff = xsumff + xff;            // calculate sigma(xi)
0618     ysumff = ysumff + yff;            // calculate sigma(yi)
0619     zsumff = zsumff + zff;            // calculate sigma(yi)
0620     
0621     x2sumff = x2sumff + xff*xff;      // calculate sigma(x^2i)
0622     y2sumff = y2sumff + yff*yff;      // calculate sigma(x^2i)
0623     z2sumff = z2sumff + zff*zff;      // calculate sigma(x^2i)
0624     
0625     xysumff = xysumff + xff * yff;    // calculate sigma(xi*yi)
0626     xzsumff = xzsumff + xff * zff;    // calculate sigma(xi*yi)
0627     yzsumff = yzsumff + yff * zff;    // calculate sigma(xi*yi)
0628     
0629     //  trkclusters.push_back(*clustertrk);
0630     //trk_nclus++;
0631     if(ptz<zfirst){zfirst = ptz;}
0632     if(ptz<zlast){zlast = ptz;}
0633       }
0634     }
0635     if(Verbosity()>0){
0636       cout << "rescheck done " <<  endl;
0637       cout << "number trk clus is " << trkclusters.size() << " " << trk_nclus<< endl;
0638     }
0639     
0640     double denominatorxff = ((x2sumff * nhitff) - (xsumff*xsumff));
0641     double slopexyff = (xysumff* nhitff - xsumff * ysumff) / denominatorxff;   // calculate slope
0642     double interceptxyff = (x2sumff * ysumff - xsumff * xysumff) / denominatorxff;  // calculate intercept
0643     
0644     double denominatorxzff = ((z2sumff * nhitff) - (zsumff*zsumff));
0645     double slopexzff = (xzsumff* nhitff - zsumff * xsumff) / denominatorxzff;   // calculate slope
0646     double interceptxzff = (z2sumff * xsumff - zsumff * xzsumff) / denominatorxzff;  // calculate intercept
0647     
0648     double denominatoryzff = ((z2sumff * nhitff) - (zsumff*zsumff));
0649     double slopeyzff = (yzsumff* nhitff - zsumff * ysumff) / denominatoryzff;   // calculate slope
0650     double interceptyzff = (z2sumff * ysumff - zsumff * yzsumff) / denominatoryzff;  // calculate intercept
0651     
0652     _ntp_trk->Fill(_nevent,tint,tsl,interceptxyff,slopexyff,trk_nclus);
0653     
0654     //  _ntp_trk_clus = new TNtuple("ntp_trk_clus", "laser clus","ev:x:y:z:adc:maxadc:size:lsize:phisize:tsize:intfit:slfit:nclus:zfirst:zlast");
0655     //_ntp_trk_hit = new TNtuple("ntp_trk_hit", "laser hitinfo","ev:x:y:z:adc:maxadc:lay:phibin:zbin:intfit:slfit:nclus:zfirst:zlast");
0656     
0657     //clean based on z position
0658     xsumff = 0;
0659     ysumff = 0;
0660     zsumff = 0;
0661     x2sumff = 0;
0662     y2sumff = 0;
0663     z2sumff = 0;
0664     xysumff = 0;
0665     xzsumff = 0;
0666     yzsumff = 0;
0667     
0668     // The shortest distance of a point from a circle is along the radial; line from the circle center to the point
0669     
0670     double axz =  slopexzff;//-slope;
0671     double bxz = -1.0;
0672     double cxz =  interceptxzff;//-intercept;
0673     
0674     double ayz =  slopeyzff;//-slope;
0675     double byz = -1.0;
0676     double cyz =  interceptyzff;//-intercept;
0677     
0678     for(vector<pointInd>::iterator clustertrk = lineclusters.begin();clustertrk!=lineclusters.end();++clustertrk){
0679       
0680       float ptx = clustertrk->first.get<0>();
0681       float pty = clustertrk->first.get<1>();
0682       TrkrDefs::cluskey tcind = clustertrk->second;
0683       LaserCluster *las_clus = _cluster_map->findCluster(tcind);      
0684       float ptz = las_clus->getZ();
0685       
0686       float res_xz =  std::abs(axz*ptz+bxz*ptx+cxz)/sqrt(axz*axz+bxz*bxz);
0687       float res_yz =  std::abs(ayz*ptz+byz*pty+cyz)/sqrt(ayz*ayz+byz*byz);
0688       //float res = std::abs(tsl*ptx+(-1*pty)+tint)/sqrt((tsl*tsl+1));
0689       if(Verbosity()>0){
0690     cout << " x: " << ptx << " z: " << ptz << " res_xz " << res_xz << endl;
0691     cout << " y: " << pty << " z: " << ptz << " res_yz " << res_yz << endl;
0692       }
0693       if(res_yz<4&&res_xz<4){
0694     float boxx = ptx;
0695     float boxy = pty;
0696     float boxz = ptz;
0697     if(Verbosity()>0){std::cout << "fprobe x:"  << boxx << " y: " << boxy <<  " z: " << boxz << std::endl;}
0698     
0699     nhitff++;
0700     float xff = boxx;
0701     float yff = boxy;
0702     float zff = boxz;
0703     
0704     xsumff = xsumff + xff;            // calculate sigma(xi)
0705     ysumff = ysumff + yff;            // calculate sigma(yi)
0706     zsumff = zsumff + zff;            // calculate sigma(yi)
0707     
0708     x2sumff = x2sumff + xff*xff;      // calculate sigma(x^2i)
0709     y2sumff = y2sumff + yff*yff;      // calculate sigma(x^2i)
0710     z2sumff = z2sumff + zff*zff;      // calculate sigma(x^2i)
0711     
0712     xysumff = xysumff + xff * yff;    // calculate sigma(xi*yi)
0713     xzsumff = xzsumff + xff * zff;    // calculate sigma(xi*yi)
0714     yzsumff = yzsumff + yff * zff;    // calculate sigma(xi*yi)
0715     
0716     trkclusters.push_back(*clustertrk);
0717     trk_nclus++;
0718     if(ptz<zfirst){zfirst = ptz;}
0719     if(ptz<zlast){zlast = ptz;}
0720       }
0721     }
0722     if(Verbosity()>0){
0723       cout << "rescheck done " <<  endl;
0724       cout << "number trk clus is " << trkclusters.size() << " " << trk_nclus<< endl;
0725     }
0726     
0727     denominatorxff = ((x2sumff * nhitff) - (xsumff*xsumff));
0728     slopexyff = (xysumff* nhitff - xsumff * ysumff) / denominatorxff;   // calculate slope
0729     interceptxyff = (x2sumff * ysumff - xsumff * xysumff) / denominatorxff;  // calculate intercept
0730     
0731     denominatorxzff = ((z2sumff * nhitff) - (zsumff*zsumff));
0732     slopexzff = (xzsumff* nhitff - xsumff * zsumff) / denominatorxzff;   // calculate slope
0733     interceptxzff = (z2sumff * xsumff - zsumff * xzsumff) / denominatorxzff;  // calculate intercept
0734     
0735     denominatoryzff = ((z2sumff * nhitff) - (zsumff*zsumff));
0736     slopeyzff = (yzsumff* nhitff - zsumff * ysumff) / denominatoryzff;   // calculate slope
0737     interceptyzff = (z2sumff * ysumff - zsumff * yzsumff) / denominatoryzff;  // calculate intercept
0738     
0739     _ntp_trk->Fill(_nevent,tint,tsl,interceptxyff,slopexyff,trk_nclus);
0740     
0741     //  _ntp_trk_clus = new TNtuple("ntp_trk_clus", "laser clus","ev:x:y:z:adc:maxadc:size:lsize:phisize:tsize:intfit:slfit:nclus:zfirst:zlast");
0742     //_ntp_trk_hit = new TNtuple("ntp_trk_hit", "laser hitinfo","ev:x:y:z:adc:maxadc:lay:phibin:zbin:intfit:slfit:nclus:zfirst:zlast"
0743     
0744     for(vector<pointInd>::iterator trkclusiter = trkclusters.begin();trkclusiter!=trkclusters.end();++trkclusiter){
0745       TrkrDefs::cluskey tcind = trkclusiter->second;
0746       LaserCluster *las_clus = _cluster_map->findCluster(tcind);
0747       float tcx = las_clus->getX();
0748       float tcy = las_clus->getY();
0749       float tcz = las_clus->getZ();
0750       
0751       unsigned int adcSum = 0.0;
0752       
0753       unsigned int maxAdc = 0.0;
0754       int iphimin = 6666, iphimax = -1;
0755       int ilaymin = 6666, ilaymax = -1;
0756       float itmin = 66666666.6, itmax = -6666666666.6;
0757       
0758       unsigned int nHits = las_clus->getNhits();
0759       if(Verbosity()>0){std::cout << "nhits " << nHits << std::endl;
0760       if(adcSum>1)
0761     puts("strange");
0762       }
0763       for(unsigned int i  = 0; i< nHits;i++){
0764     float coords[3] = {las_clus->getHitX(i), las_clus->getHitY(i), las_clus->getHitZ(i)};
0765     int iphi = las_clus->getHitIPhi(i);
0766     if(Verbosity()>0){std::cout << " iphi " << iphi << std::endl;
0767     }
0768     int it = las_clus->getHitIT(i);
0769     int ilay = las_clus->getHitLayer(i);
0770     unsigned int adc = las_clus->getHitAdc(i);
0771     if(iphi<iphimin){iphimin = iphi;}
0772     if(iphi>iphimax){iphimax = iphi;}
0773     if(ilay<ilaymin){ilaymin = ilay;}
0774     if(ilay>ilaymax){ilaymax = ilay;}
0775     if(it<itmin){itmin = it;}
0776     if(it>itmax){itmax = it;}
0777     adcSum += adc;
0778     if(adc>maxAdc){
0779       maxAdc = adc;
0780     }
0781     //  printf("test %d %f %f %f %d %d %d %d %f %f %d %f %f\n", _nevent, coords[0],coords[1],coords[2],adc,ilay,iphi,it,interceptxyff,slopexyff,trk_nclus,zfirst,zlast);
0782     
0783     m_nevent = _nevent;
0784     m_hitx = coords[0];
0785     m_hity = coords[1];
0786     m_hitz = coords[2];
0787     m_hitadc = adc;
0788     m_hitlayer = ilay;
0789     m_hitpad = iphi;
0790     m_hittbin = it;
0791     m_hittbin = process_side;
0792     m_interxy = interceptxyff;
0793     m_slopexy = slopexyff;
0794     m_interxz = interceptxzff;
0795     m_slopexz = slopexzff;
0796     m_interyz = interceptyzff;
0797     m_slopeyz = slopeyzff;
0798     m_nclus = trk_nclus;
0799     m_zfirst = zfirst;
0800     m_zlast = zlast;
0801     
0802     m_hittree->Fill(); 
0803       } 
0804       
0805       int phisize =  iphimax - iphimin + 1;
0806       int lsize =  ilaymax - ilaymin + 1;
0807       int tsize = itmax - itmin +1;
0808       
0809       m_clux = tcx;
0810       m_cluy = tcy;
0811       m_cluz = tcz;
0812 
0813       m_cluadc = adcSum;
0814       m_clumaxadc = maxAdc;
0815       m_size = nHits;
0816       m_sizel = lsize;
0817       m_sizephi =  phisize;
0818       m_sizet = tsize;
0819       
0820       m_clustree->Fill();
0821     }
0822   }
0823   }
0824   if(Verbosity()>0){ cout << "number of seeds is " << numberofseeds << endl;}
0825   if(_write_ntp){
0826     _ntp_max->Fill(_nevent,intmin,intmax,slmin,slmax);
0827   }  
0828  
0829   if(Verbosity()>0){cout << "return " << numberofseeds << endl;}
0830   return Fun4AllReturnCodes::EVENT_OK;
0831 }
0832 
0833 
0834 int PHLineLaserReco::Setup(PHCompositeNode* topNode)
0835 {
0836   if(Verbosity()>0){
0837     cout << "Called Setup" << endl;
0838     cout << "topNode:" << topNode << endl;
0839   }
0840   // PHTrackSeeding::Setup(topNode);
0841   //  int ret = GetNodes(topNode);
0842   //return ret;
0843   GetNodes(topNode);
0844   return Fun4AllReturnCodes::EVENT_OK;
0845 }
0846 
0847 int PHLineLaserReco::End(PHCompositeNode* topNode)
0848 {
0849   if(topNode==nullptr){
0850     std::cout << PHWHERE << "No topNode, bailing."
0851           << std::endl;
0852     return Fun4AllReturnCodes::ABORTEVENT;
0853   }
0854     
0855   if(_write_ntp){
0856     _tfile->cd();
0857     _ntp_cos->Write();
0858     _ntp_stub->Write();
0859     _ntp_max->Write();
0860     _ntp_trk->Write();
0861     m_hittree->Write();
0862     m_clustree->Write();
0863     _tfile->Close();
0864     if(Verbosity()>0){ cout << "Called End " << endl;}
0865   }
0866   return Fun4AllReturnCodes::EVENT_OK;
0867 }
0868 
0869
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!