Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 /*!
0002  *  \file PHStreakFinder.cc
0003  
0004  *  \author Christof Roland 
0005  */
0006 
0007 //begin
0008 
0009 //#ifndef PHCOSMICSFILTER_H
0010 //#define PHCOSMICSFILTER_H
0011 
0012 #include <fun4all/SubsysReco.h>
0013 #include <trackbase/TrkrCluster.h>
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase_historic/TrackSeed.h>
0016 #include <trackbase_historic/TrackSeedContainer.h>
0017 #include <trackbase_historic/TrackSeedContainer_v1.h>
0018 #include <trackbase_historic/TrackSeed_v2.h>
0019 
0020 // Helix Hough + Eigen includes (hidden from rootcint)
0021 #if !defined(__CINT__) || defined(__CLING__)
0022 
0023 #include <Eigen/Core>                              // for Matrix
0024 #endif
0025 
0026 
0027 
0028 #if !defined(__CINT__) || defined(__CLING__)
0029 //BOOST for combi seeding
0030 #include <boost/geometry.hpp>
0031 #include <boost/geometry/geometries/box.hpp>
0032 #include <boost/geometry/geometries/point.hpp>
0033 
0034 #include <boost/geometry/index/rtree.hpp>
0035 #endif
0036 
0037 // standard includes
0038 #include <cfloat>
0039 #include <iostream>                                // for operator<<, basic_...
0040 #include <list>
0041 #include <map>
0042 #include <memory>
0043 #include <set>                                     // for set
0044 #include <string>                                  // for string
0045 #include <utility>                                 // for pair
0046 #include <vector>
0047 
0048 // forward declarations
0049 class BbcVertexMap;
0050 
0051 class PHCompositeNode;
0052 
0053 class PHG4CellContainer;
0054 class PHG4CylinderGeomContainer;
0055 class PHG4Particle;
0056 class PHG4TruthInfoContainer;
0057 class PHTimer;
0058 
0059 class sPHENIXSeedFinder;
0060 
0061 class SvtxClusterMap;
0062 class SvtxCluster;
0063 class SvtxTrackMap;
0064 class SvtxTrack;
0065 class SvtxTrackState;
0066 class SvtxVertexMap;
0067 class SvtxHitMap;
0068 
0069 class TNtuple;
0070 class TFile;
0071 
0072 #include "PHStreakFinder.h"
0073 
0074 #include "AssocInfoContainer.h"                         // for AssocInfoCont...
0075 
0076 // trackbase_historic includes
0077 #include <trackbase/TrackFitUtils.h>
0078 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0079 #include <trackbase/TrkrClusterContainer.h>
0080 #include <trackbase/TrkrDefs.h>  // for getLayer, clu...
0081 #include <trackbase/TrkrClusterHitAssoc.h>
0082 #include <trackbase/TrkrClusterIterationMapv1.h>
0083 #include <trackbase_historic/TrackSeedContainer.h>
0084 #include <trackbase_historic/TrackSeed_v2.h>
0085 #include <trackbase/LaserCluster.h>
0086 #include <trackbase/LaserClusterContainer.h>
0087 #include <trackbase/LaserClusterContainerv1.h>
0088 #include <trackbase/LaserClusterv1.h>
0089 
0090 // sPHENIX Geant4 includes
0091 #include <g4detectors/PHG4CylinderCellGeom.h>
0092 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0093 #include <g4detectors/PHG4CylinderGeom.h>
0094 #include <g4detectors/PHG4CylinderGeomContainer.h>
0095 
0096 // sPHENIX includes
0097 #include <fun4all/Fun4AllReturnCodes.h>
0098 
0099 #include <phool/PHTimer.h>                              // for PHTimer
0100 #include <phool/getClass.h>
0101 
0102 #include <Eigen/Core>                  // for Matrix
0103 #include <Eigen/Dense>
0104 
0105 
0106 
0107 //ROOT includes for debugging
0108 #include <TFile.h>
0109 #include <TNtuple.h>
0110 #include <TAxis.h>
0111 #include <TGraph.h>
0112 //BOOST for combi seeding
0113 #include <boost/geometry/geometries/box.hpp>
0114 #include <boost/geometry/geometries/point.hpp>
0115 #include <boost/geometry/index/rtree.hpp>
0116 
0117 using point = bg::model::point<float, 3, bg::cs::cartesian>;
0118 //typedef bg::model::point<float, 3, bg::cs::cartesian> point;
0119 using box = bg::model::box<point>;
0120 //typedef bg::model::box<point> box;
0121 using pointKey = std::pair<point, int>;
0122 using stubKey = std::pair<point, float>;
0123 //typedef std::pair<point, TrkrDefs::cluskey> pointKey;
0124 
0125 using pointrtree = bgi::rtree<pointKey, bgi::quadratic<16>>;
0126 using stubrtree = bgi::rtree<pointKey, bgi::quadratic<16>>;
0127 
0128 // standard includes
0129 #include <TH1.h>
0130 //#include <stdio.h>      /* printf */
0131 //#include <cmath.h>       /* copysign */
0132 #include <algorithm>
0133 #include <climits>                                     // for UINT_MAX
0134 #include <cmath>
0135 #include <iostream>
0136 #include <memory>
0137 #include <set>                                          // for set
0138 #include <tuple>
0139 #include <utility>   
0140                                    // for pair, make_pair
0141 
0142 #define LogDebug(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0143 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
0144 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
0145 
0146 
0147 using namespace std;
0148 //using namespace ROOT::Minuit2;
0149 //namespace bg = boost::geometry;
0150 //namespace bgi = boost::geometry::index;
0151 
0152 //typedef uint64_t cluskey;
0153 
0154 //vector<LaserCluster*> clusterpoints;
0155 using myvec = vector<float>;
0156 
0157 PHStreakFinder::PHStreakFinder(const std::string& name, const string& filename)
0158   : SubsysReco(name)
0159 , _filename(filename)
0160 {
0161   
0162 }
0163 
0164 int PHStreakFinder::GetNodes(PHCompositeNode* topNode)
0165 {
0166   tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0167   if (!tGeometry)
0168   {
0169     std::cout << PHWHERE << "No acts reco geometry, bailing."
0170               << std::endl;
0171     return Fun4AllReturnCodes::ABORTEVENT;
0172   }
0173   _cluster_map = findNode::getClass<LaserClusterContainerv1>(topNode, "LASER_CLUSTER");
0174   if (!_cluster_map)
0175   {
0176     std::cout << PHWHERE << "No cluster container on node tree, bailing."
0177               << std::endl;
0178     return Fun4AllReturnCodes::ABORTEVENT;
0179   }
0180   return Fun4AllReturnCodes::EVENT_OK;
0181 }
0182 
0183 //____________________________________________________________________________..
0184 int PHStreakFinder::Init(PHCompositeNode* topNode)
0185 {
0186   if(topNode==nullptr){
0187     std::cout << PHWHERE << "No topNode, bailing."
0188           << std::endl;
0189     return Fun4AllReturnCodes::ABORTEVENT;
0190   }
0191   _nevent = 0;
0192     _tfile = new TFile(_filename.c_str(), "RECREATE");
0193   _ntp_ev = new TNtuple("ntp_ev", "cos event info","ev:ev2");
0194   _ntp_clu = new TNtuple("ntp_clu", "cos event info","ev:x:y:z");
0195   _ntp_clutrk = new TNtuple("ntp_clutrk", "cos event info","ev:ntrk:x:y:z:nclus:nclus0:nclus1:r:phi:dedx:intxz:intyz:intrz:slxz:slyz:slrz:mresx:mresy:zmin:zmax");
0196   _ntp_stub = new TNtuple("ntp_stub", "cos stub info","ev:intxz:slxz:intyz:slyz:tanxz:tanyz:num");
0197   _ntp_trk = new TNtuple("ntp_trk", "cos stub info","ev:ntrk:nclus:nclus0:nclus1:r:phi:dedx:intxz:intyz:intrz:slxz:slyz:slrz:mresx:mresy:zmin:zmax");
0198 
0199   return Fun4AllReturnCodes::EVENT_OK;
0200 }
0201 
0202 //____________________________________________________________________________..
0203 int PHStreakFinder::InitRun(PHCompositeNode* topNode)
0204 {
0205   int ret = GetNodes(topNode);
0206   if (ret != Fun4AllReturnCodes::EVENT_OK)
0207   {
0208     return ret;
0209   }
0210   ret = createNodes(topNode);
0211 
0212 
0213   return ret;
0214 }
0215 int PHStreakFinder::createNodes(PHCompositeNode* topNode)
0216 {
0217   PHNodeIterator iter(topNode);
0218 
0219   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0220 
0221   if (!dstNode)
0222   {
0223     std::cerr << "DST node is missing, quitting" << std::endl;
0224     throw std::runtime_error("Failed to find DST node in PHStreakFinder::createNodes");
0225   }
0226 
0227   PHNodeIterator dstIter(dstNode);
0228   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0229 
0230   if (!svtxNode)
0231   {
0232     svtxNode = new PHCompositeNode("SVTX");
0233     dstNode->addNode(svtxNode);
0234   }
0235 
0236   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
0237   if (!m_seedContainer)
0238   {
0239     m_seedContainer = new TrackSeedContainer_v1;
0240     PHIODataNode<PHObject>* trackNode =
0241         new PHIODataNode<PHObject>(m_seedContainer, m_trackMapName, "PHObject");
0242     svtxNode->addNode(trackNode);
0243   }
0244   if (m_seedContainer){
0245     cout << "SEED CONTAINER CREATED" << endl;
0246   }
0247   return Fun4AllReturnCodes::EVENT_OK;
0248 }
0249 
0250 double PHStreakFinder::phiadd(double phi1, double phi2){
0251   double s = phi1+phi2;
0252   if(s>2*M_PI) {return s-2*M_PI;}
0253   else if(s<0) {return s+2*M_PI;}
0254   else {return s;}
0255 }
0256 
0257 
0258 double PHStreakFinder::phidiff(double phi1, double phi2){
0259   double d = phi1-phi2;
0260   if(d>M_PI) {return d-2*M_PI;}
0261   else if(d<-M_PI) {return d+2*M_PI;}
0262   else {return d;}
0263 }
0264 
0265 void PHStreakFinder::get_stub(const pointrtree &search_rtree, float pointx, float pointy, float pointz, int &count, double &xzslope, double &xzintercept, double &yzslope, double &yzintercept){ //NOLINT
0266   float m1_dx = 2;
0267   float m1_dy = 2;
0268   float m1_dz = 8;
0269   vector<pointKey> boxclusters;
0270   search_rtree.query(bgi::intersects(box(point(pointx-m1_dx,pointy-m1_dy,pointz-m1_dz),
0271                   point(pointx+m1_dx,pointy+m1_dy,pointz+m1_dz))),std::back_inserter(boxclusters));
0272   
0273   int nbox = boxclusters.size();
0274   
0275   int nhit = 0;
0276   double xsum = 0;
0277   double x2sum = 0;
0278   double ysum = 0;
0279   double y2sum = 0;
0280   double z2sum = 0;
0281 
0282   double zsum = 0;
0283   double xzsum = 0;
0284   double yzsum = 0;
0285 
0286   if(nbox>=2){
0287     for(vector<pointKey>::iterator pbox = boxclusters.begin();pbox!=boxclusters.end();++pbox){
0288       float boxx = pbox->first.get<0>();
0289       float boxy = pbox->first.get<1>();
0290       float boxz = pbox->first.get<2>();
0291       //      std::cout << "fprobe: " << count << " x:"  << boxx << " y: " << boxy << std::endl;
0292       nhit++;
0293       float x = boxx;
0294       float y = boxy;
0295       float z = boxz;
0296 
0297       xsum = xsum + x;            // calculate sigma(xi)
0298       ysum = ysum + y;            // calculate sigma(xi)
0299       zsum = zsum + z;            // calculate sigma(yi)
0300 
0301       x2sum = x2sum + x*x;  // calculate sigma(x^2i)
0302       y2sum = y2sum + y*y;  // calculate sigma(x^2i)
0303       z2sum = z2sum + z*z;  // calculate sigma(x^2i)
0304 
0305       xzsum = xzsum + x * z;      // calculate sigma(xi*yi)
0306       yzsum = yzsum + y * z;      // calculate sigma(xi*yi)
0307     }
0308   }
0309   
0310   const double denominator = ((z2sum * nhit) - (zsum*zsum));
0311   //  const double denominator_yz = ((y2sum * nhit) - (ysum*ysum));
0312 
0313   xzslope = (xzsum * nhit - zsum * xsum) / denominator;   // calculate slope
0314   xzintercept = (z2sum * xsum - zsum * xzsum) / denominator;  // calculate intercept
0315 
0316   yzslope = (yzsum * nhit - zsum * ysum) / denominator;   // calculate slope
0317   yzintercept = (z2sum * ysum - zsum * yzsum) / denominator;  // calculate intercept
0318 
0319   count = nhit;
0320 }
0321 
0322 int PHStreakFinder::process_event(PHCompositeNode* topNode)
0323 {
0324   if(topNode==nullptr){
0325     std::cout << PHWHERE << "No topNode, bailing."
0326           << std::endl;
0327     return Fun4AllReturnCodes::ABORTEVENT;
0328   }
0329     
0330   _nevent++;
0331   _ntp_ev->Fill(_nevent,_nevent);
0332   int ntrack = 0;
0333   std::vector<TrackSeed_v2> clean_chains;
0334   //Fill rtree
0335   bgi::rtree<pointKey, bgi::quadratic<16> > rtree;
0336 
0337   auto range = _cluster_map->getClusters();
0338   for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
0339     LaserCluster *cluster = clusIter->second;
0340     _ntp_clu->Fill(_nevent,cluster->getX() , cluster->getY(), cluster->getZ());
0341     rtree.insert(std::make_pair(point(cluster->getX() , cluster->getY(), cluster->getZ()),std::distance(range.first,clusIter) ));
0342       
0343   }
0344 
0345   //Get all clusters from rtree, fit stubs around clusters, fill stub tree
0346   // search horizontal strips in x - y search grid
0347 
0348   for(float ix = -80;ix<80;ix+=4){
0349     for(float iy = -80;iy<80;iy+=4){ 
0350       float ir = sqrt(ix*ix+iy*iy);
0351       if(ir<30||ir>78){
0352     continue;
0353       }
0354       vector<pointKey> allclusters;
0355       
0356       rtree.query(bgi::intersects(box(point(ix-2,iy-2,-110),point(ix+2,iy+2,110))),std::back_inserter(allclusters));
0357       //if(Verbosity()>0)
0358       if(allclusters.size()>20){
0359     //cout << " ev: " << _nevent << " ix = " << ix << " iy: " << iy << " number clus is " << allclusters.size() << endl;
0360       
0361     
0362     bgi::rtree<stubKey, bgi::quadratic<16> > xzrtree_stub;
0363     bgi::rtree<stubKey, bgi::quadratic<16> > yzrtree_stub;
0364     
0365     for(vector<pointKey>::iterator cluster = allclusters.begin();cluster!=allclusters.end();++cluster){
0366       float pointx = cluster->first.get<0>();
0367       float pointy = cluster->first.get<1>();
0368       float pointz = cluster->first.get<2>();
0369       int fcount = 0;
0370       double xzslope = 0;
0371       double xzintercept = 0;
0372       double yzslope = 0;
0373       double yzintercept = 0;
0374       //calc slope and intersect from 5 cluster stubs
0375       get_stub(rtree, pointx, pointy, pointz, fcount, xzslope, xzintercept, yzslope, yzintercept);
0376       //int num = 0;
0377       if(finite(xzslope)&&finite(yzslope)){
0378         xzrtree_stub.insert(std::make_pair(point(xzintercept , xzslope, yzintercept), yzslope));
0379         //  cout << " inserting: " << "xzint: " << xzintercept << "xzsl: " << xzslope << " yzint: " << yzintercept << " yzsl: " << yzslope <<endl;
0380         float tana_xz = atan(xzslope);
0381         float tana_yz = atan(yzslope);
0382         //ev:intxz:slxz:intyz:slyz:tanxz:tanyz:xo
0383         if(_write_ntp)_ntp_stub->Fill(_nevent,xzintercept,xzslope,yzintercept,yzslope,tana_xz,tana_yz,fcount);
0384         
0385       }   
0386     }
0387     //find good stub candidates
0388     map <int,pair<float,float>> xzouttrkmap;
0389     map <int,myvec> xztrkmap;
0390     //  int nout_xz = 0;
0391     float intmax = 80;
0392     float intmin = -80;
0393     float slmax = 0.2;
0394     float slmin = -0.2;
0395     float int_width = (intmax - intmin)/20;
0396     float sl_width = (slmax - slmin)/8;
0397     cout << "number xzstub is " << xzrtree_stub.size() << endl;
0398     
0399     std::vector<stubKey> xzallstubs;
0400     xzrtree_stub.query(bgi::intersects(box(point(intmin,slmin,intmin),point(intmax,slmax,intmax))),std::back_inserter(xzallstubs));
0401     if(xzallstubs.size()==0){cout<<"nostubs" <<endl; break;}
0402 
0403 
0404     for(vector<stubKey>::iterator stub = xzallstubs.begin();stub!=xzallstubs.end();++stub){
0405       float pxzint = stub->first.get<0>();
0406       float pxzsl = stub->first.get<1>();
0407       float pyzint = stub->first.get<2>();
0408       float pyzsl =  stub->second;
0409       //cout << "xzint: " << pxzint << "xzsl: " << pxzsl << " yzint: " << pyzint << " yzsl: " << pyzsl <<endl;
0410       if(abs(pyzsl)<0.2){
0411         vector<stubKey> trkcand;
0412         xzrtree_stub.query(bgi::intersects(box(point(pxzint-int_width,pxzsl-sl_width,pyzint-int_width),point(pxzint+int_width,pxzsl+sl_width,pyzint+int_width))),std::back_inserter(trkcand));
0413       
0414         int ntrk = trkcand.size();
0415         int count = 0;
0416         float xzintsum = 0;
0417         float xzslsum = 0;
0418         float yzintsum = 0;
0419         float yzslsum = 0;
0420         if(ntrk>=5){
0421           for(vector<stubKey>::iterator ptrk = trkcand.begin();ptrk!=trkcand.end();++ptrk){
0422         float trkxzint = ptrk->first.get<0>();
0423         float trkxzsl = ptrk->first.get<1>();
0424         float trkyzint = ptrk->first.get<2>();
0425         float trkyzsl = ptrk->second;
0426         //    if(Verbosity()>0)
0427         /*cout<< " ev: " << _nevent
0428             << " stub " << ntrk
0429             << " xzint: " << trkxzint 
0430             << " xzsl: "  << trkxzsl 
0431             << " yzint: " << trkyzint 
0432             << " yzsl: "  << trkyzsl 
0433             << endl;
0434         */
0435         xzintsum = xzintsum + trkxzint;            // calculate sigma(xi)
0436         xzslsum = xzslsum + trkxzsl;            // calculate sigma(yi)
0437         yzintsum = yzintsum + trkyzint;            // calculate sigma(xi)
0438         yzslsum = yzslsum + trkyzsl;            // calculate sigma(yi)
0439         count++;
0440           }
0441           float mxzint = (xzintsum/count);
0442           float mxzsl  = (xzslsum/count);
0443           float myzint = (yzintsum/count);
0444           float myzsl  = (yzslsum/count);
0445           myvec outarr = {mxzint,mxzsl,myzint,myzsl} ;
0446           /*          outarr.push_back(mxzint);
0447           outarr.push_back(mxzsl);
0448           outarr.push_back(myzint);
0449           outarr.push_back(myzsl);
0450           */
0451           xztrkmap[ntrk] = outarr;
0452           //.insert(std::pair<int,vector<float>(ntrk,outarr)); 
0453           
0454           if(Verbosity()>0) {
0455           cout<< " sxz tub in box " << ntrk
0456           << " xzint: " << pxzint 
0457           << " xzsl: " << pxzsl 
0458           << " yzint: " << pyzint 
0459           << " yzsl: " << pyzsl 
0460           << " intwi: " << int_width 
0461           << " slwi: " << sl_width
0462           << " mxzint " << mxzint
0463           << " mxzsl " << mxzsl
0464           << " myzint " << myzint
0465           << " myzsl " << myzsl
0466           << endl;
0467           }
0468         }
0469       }
0470     }
0471     float meanxzint = 0;
0472     float meanxzsl = 0;
0473     float meanyzint = 0;
0474     float meanyzsl = 0;
0475 //  int meancnt = 0;
0476     for(auto ptrk = xztrkmap.begin();ptrk!=xztrkmap.end();++ptrk){
0477       int nstub = ptrk->first;
0478       myvec value = ptrk->second;
0479 //    meancnt++;
0480       meanxzint = value[0];
0481       meanxzsl  = value[1];
0482       meanyzint = value[2];
0483       meanyzsl  = value[3];
0484       //      for (const auto& [nstub, value[4]] : xztrkmap)
0485        if(Verbosity()>0) {
0486          std::cout <<" nstub: " << nstub
0487             << "xzint " << value[0] 
0488             << "xzsl  " << value[1] 
0489             << "yzint " << value[2] 
0490             << "yzsl " << value[3] 
0491             << std::endl;
0492        }
0493     }
0494     //collect hits and fill ntuple
0495     
0496     vector<pointKey> trkclusters;
0497     vector<pointKey> lineclusters;
0498     rtree.query(bgi::intersects(box(point(-80,-80,-110),point(80,80,110))),std::back_inserter(lineclusters));
0499     // if(Verbosity()>0) 
0500      if(Verbosity()>0) {
0501        cout << "number line clus is " << lineclusters.size() << endl;
0502      }
0503     float meanresx = 0;
0504     float meanresy = 0;
0505       
0506     for(vector<pointKey>::iterator clustertrk = lineclusters.begin();clustertrk!=lineclusters.end();++clustertrk){
0507       
0508       float ptx = clustertrk->first.get<0>();
0509       float pty = clustertrk->first.get<1>();
0510       float ptz = clustertrk->first.get<2>();
0511       
0512       float refx = meanxzint + ptz*meanxzsl;
0513       float refy = meanyzint + ptz*meanyzsl;
0514       float resx =  abs(ptx-refx);
0515       float resy =  abs(pty-refy);
0516       //      if(Verbosity()>0)
0517       if((resx<2)&&(resy<2)){
0518         meanresx += resx;
0519         meanresy += resy;
0520         if(Verbosity()>0) { cout << " x: " << ptx << " | " << refx << " y: " << pty << " | " << refy << " resx " << resx << " resy:" << resy << endl;}
0521         trkclusters.push_back(*clustertrk);
0522       }
0523     }
0524     meanresx/= trkclusters.size();
0525     meanresy/= trkclusters.size();
0526     int nclus = 0;
0527     int nclus0 = 0;
0528     int nclus1 = 0;
0529     float zmin = 10000;
0530     float zmax = -10000;
0531     std::vector<float> dedxlist;
0532     
0533     //cout << "number trk clus is " << trkclusters.size() << endl;
0534     if(trkclusters.size()>=_min_nclusters){
0535       //cout << "setting keep true: " << trkclusters.size() << " > " << _min_nclusters << endl;
0536       for (const auto& trkclu : trkclusters){
0537         //      std::cout <<" yyz cand out ev: " << _nevent << '[' << xzkey << "] = " << yzvalue.first << " | " << xzvalue.second << std::endl;
0538         // for(vector<pointKey>::iterator clustertrk = trkclusters.begin();clustertrk!=trkclusters.end();++clustertrk){
0539         float tptx = trkclu.first.get<0>();
0540         float tpty = trkclu.first.get<1>();
0541         float tptz = trkclu.first.get<2>();
0542         //      float ptx = clustertrk->first.get<0>();
0543         //float pty = clustertrk->first.get<1>();
0544         // float ptz = clustertrk->first.get<2>();
0545         nclus++;
0546         if(tptz<0){
0547           nclus0++;
0548         }
0549         if(tptz>0){
0550           nclus1++;
0551         }
0552         if(tptz<zmin){
0553           zmin = tptz;
0554         }
0555         if(tptz>zmax){
0556           zmax = tptz;
0557         }
0558         float adc = trkclu.second;
0559         
0560         float thick = 1.098;
0561         float r  = sqrt(tptx*tptx+tpty*tpty);
0562         if(r>30&&r<=40.5){
0563           thick = 0.57;
0564         }else if(r>40.5&&r<60){
0565           thick = 1.02;
0566         } 
0567         dedxlist.push_back(adc/thick);
0568         //calc dedx
0569         //_ntp_clutrk->Fill(_nevent,tptx , tpty, tptz);
0570       }
0571       sort(dedxlist.begin(), dedxlist.end());
0572       int trunc_min = 0;
0573       int trunc_max = (int)dedxlist.size()*0.7;
0574       float sumdedx = 0;
0575       int ndedx = 0;
0576       for(int j = trunc_min; j<=trunc_max;j++){
0577         sumdedx+=dedxlist.at(j);
0578         ndedx++;
0579       }
0580       sumdedx/=ndedx;
0581       TVector2 vec(meanxzint,meanyzint);
0582       // ev:nclus:nclus0:nclus1:r:phi:dedx:intxz:intyz:intrz:slxz:slyz:slrz
0583       float r = sqrt(meanxzint*meanxzint+ meanyzint*meanyzint);
0584       float xp = meanxzint+meanxzsl*100;
0585       float xn = meanxzint+meanxzsl*-100;
0586       float yp = meanyzint+meanyzsl*100;
0587       float yn = meanyzint+meanyzsl*-100;
0588       float rp = sqrt(xp*xp+yp*yp);
0589       float rn = sqrt(xn*xn+yn*yn);
0590       float phi = vec.Phi();
0591       float dr = rp-rn;
0592       float ftX[30] = {0};
0593       int nft = 0;
0594       ftX[nft++] = _nevent;
0595       ftX[nft++] =ntrack;
0596       ftX[nft++] =nclus;
0597       ftX[nft++] =nclus0;
0598       ftX[nft++] =nclus1;
0599       ftX[nft++] =r;
0600       ftX[nft++] =phi;
0601       ftX[nft++] =sumdedx;
0602       ftX[nft++] =meanxzint;
0603       ftX[nft++] =meanyzint;
0604       ftX[nft++] =r;
0605       ftX[nft++] =meanxzsl;
0606       ftX[nft++] =meanyzsl;
0607       ftX[nft++] =dr/200;
0608       ftX[nft++] =meanresx;
0609       ftX[nft++] =meanresy;
0610       ftX[nft++] =zmin;
0611       ftX[nft++] =zmax;
0612       _ntp_trk->Fill(ftX);
0613 
0614       for (const auto& trkclu : trkclusters){
0615         float tptx = trkclu.first.get<0>();
0616         float tpty = trkclu.first.get<1>();
0617         float tptz = trkclu.first.get<2>();
0618         float fcX[30] = {0};
0619         int nfc = 0;
0620         fcX[nfc++] = _nevent;
0621         fcX[nfc++] =ntrack;
0622         fcX[nfc++] =tptx;
0623         fcX[nfc++] =tpty;
0624         fcX[nfc++] =tptz;
0625         fcX[nfc++] =nclus;
0626         fcX[nfc++] =nclus0;
0627         fcX[nfc++] =nclus1;
0628         fcX[nfc++] =r;
0629         fcX[nfc++] =phi;
0630         fcX[nfc++] =sumdedx;
0631         fcX[nfc++] =meanxzint;
0632         fcX[nfc++] =meanyzint;
0633         fcX[nfc++] =r;
0634         fcX[nfc++] =meanxzsl;
0635         fcX[nfc++] =meanyzsl;
0636         fcX[nfc++] =dr/200;
0637         fcX[nfc++] =meanresx;
0638         fcX[nfc++] =meanresy;
0639         fcX[nfc++] =zmin;
0640         fcX[nfc++] =zmax;
0641         
0642         _ntp_clutrk->Fill(fcX);
0643       }
0644       ntrack++;
0645     }
0646       }
0647     }
0648   }
0649 
0650 
0651  #ifdef KAKAKA
0652 
0653   
0654   //  int out = 0;
0655   // find clusters of slope intectept pairs
0656   map <int,pair<float,float>> xzouttrkmap;
0657   map <int,pair<float,float>> yzouttrkmap;
0658   map <int,pair<float,float>> xztrkmap;
0659   map <int,pair<float,float>> yztrkmap;
0660   int nout_xz = 0;
0661   int nout_yz = 0;
0662   float intmax = 80;
0663   float intmin = -80;
0664   float slmax = 0.5;
0665   float slmin = -0.5;
0666   float int_width = (intmax - intmin)/20;
0667   float sl_width = (slmax - slmin)/20;
0668   cout << "number xzstuc is " << xzrtree_stub.size() << endl;
0669   cout << "number yzstub is " << yzrtree_stub.size() << endl;
0670   while(xzrtree_stub.size()>10){
0671     std::vector<pointKey> xzallstubs;
0672     xzrtree_stub.query(bgi::intersects(box(point(intmin,slmin,-1),point(intmax,slmax,1))),std::back_inserter(xzallstubs));
0673     if(xzallstubs.size()==0){cout<<"nostubs" <<endl; break;}
0674 
0675 
0676     for(vector<pointKey>::iterator stub = xzallstubs.begin();stub!=xzallstubs.end();++stub){
0677       float pint = stub->first.get<0>();
0678       float psl = stub->first.get<1>();
0679       cout << "int: " << pint << "sl: " << psl << endl;
0680       vector<pointKey> trkcand;
0681       
0682       xzrtree_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));
0683       
0684       int ntrk = trkcand.size();
0685       int count = 0;
0686       float intsum = 0;
0687       float slsum = 0;
0688       if(ntrk>=5){
0689     for(vector<pointKey>::iterator ptrk = trkcand.begin();ptrk!=trkcand.end();++ptrk){
0690       float trkint = ptrk->first.get<0>();
0691       float trksl = ptrk->first.get<1>();
0692       //      if(Verbosity()>0)
0693         cout<< " ev: " << _nevent
0694         << " xz stub " << ntrk
0695         << " int: " << trkint 
0696         << " sl: "  << trksl 
0697         << endl;
0698       
0699       intsum = intsum + trkint;            // calculate sigma(xi)
0700       slsum = slsum + trksl;            // calculate sigma(yi)
0701       count++;
0702     }
0703     float mint = (intsum/count);
0704     float msl  = (slsum/count);
0705     xztrkmap[ntrk] = std::make_pair(mint,msl); 
0706       
0707     //  if(Verbosity()>0) 
0708     cout<< " sxz tub in box " << ntrk
0709         << " int: " << pint 
0710         << " sl: " << psl 
0711         << " intwi: " << int_width 
0712         << " slwi: " << sl_width
0713         << " mint " << mint
0714         << " msl " << msl
0715         << endl;
0716       }
0717     }
0718   while(yzrtree_stub.size()>10){
0719     std::vector<pointKey> yzallstubs;
0720     yzrtree_stub.query(bgi::intersects(box(point(intmin,slmin,-1),point(intmax,slmax,1))),std::back_inserter(yzallstubs));
0721     if(yzallstubs.size()==0){break;}
0722 
0723     for(vector<pointKey>::iterator stub = yzallstubs.begin();stub!=yzallstubs.end();++stub){
0724       float pint = stub->first.get<0>();
0725       float psl = stub->first.get<1>();
0726       vector<pointKey> trkcand;
0727       
0728       yzrtree_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));
0729       
0730       int ntrk = trkcand.size();
0731       int count = 0;
0732       float intsum = 0;
0733       float slsum = 0;
0734       if(ntrk>=5){
0735     for(vector<pointKey>::iterator ptrk = trkcand.begin();ptrk!=trkcand.end();++ptrk){
0736       float trkint = ptrk->first.get<0>();
0737       float trksl = ptrk->first.get<1>();
0738       //      if(Verbosity()>0)
0739         cout<< " ev: " << _nevent
0740         << " yz stub " << ntrk
0741         << " int: " << trkint 
0742         << " sl: "  << trksl 
0743         << endl;
0744       
0745       intsum = intsum + trkint;            // calculate sigma(xi)
0746       slsum = slsum + trksl;            // calculate sigma(yi)
0747       count++;
0748     }
0749     float mint = (intsum/count);
0750     float msl  = (slsum/count);
0751     yztrkmap[ntrk] = std::make_pair(mint,msl); 
0752       
0753     //  if(Verbosity()>0) 
0754     cout<< " syz tub in box " << ntrk
0755         << " int: " << pint 
0756         << " sl: " << psl 
0757         << " intwi: " << int_width 
0758         << " slwi: " << sl_width
0759         << " mint " << mint
0760         << " msl " << msl
0761         << endl;
0762       }
0763     }
0764     
0765     if( xztrkmap.size()>0){
0766       int size_before = xzrtree_stub.size();
0767       //      if(Verbosity()>0)
0768       cout << "mapend: " << xztrkmap.rbegin()->first << " | " << (xztrkmap.rbegin()->second).first << " | "  << xztrkmap.rbegin()->second.second << endl;
0769       // for (const auto& [key, value] : trkmap)
0770       //    std::cout <<" ev: " << _nevent << '[' << key << "] = " << value.first << " | " << value.second << std::endl;
0771       //remove stubs in box
0772       float rmint = xztrkmap.rbegin()->second.first;
0773       float rmsl = xztrkmap.rbegin()->second.second;
0774       vector<pointKey> rmcand;
0775       xzrtree_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));
0776       for(vector<pointKey>::iterator rmstub = rmcand.begin();rmstub!=rmcand.end();++rmstub){
0777     float rmpint = rmstub->first.get<0>();
0778     float rmpsl = rmstub->first.get<1>();
0779     if(Verbosity()>0) cout<< "    rm " <<  " int: " << rmpint 
0780                   << " sl: " << rmpsl 
0781                   << endl;
0782     
0783     //auto rmpoint = std::make_pair(point(rmpint , rmpsl, 0.0), 0);
0784     //rtree_stub.remove(rmpoint);
0785     //  pointKey match;
0786     xzrtree_stub.remove(*rmstub);
0787     /*  if (!rtree_stub.query(bgi::nearest(rmpoint, 1), &match))
0788         continue;
0789     */
0790     /*
0791       if (bg::distance(rmpoint, match.first) <= 0.5) {
0792       rtree_stub.remove(match);
0793       
0794     */
0795       }
0796       int size_after = xzrtree_stub.size();
0797       if(size_before == size_after) {break;}
0798       xzouttrkmap[nout_xz++] = std::make_pair(rmint,rmsl);
0799       cout << " tree sixe after remove: " << xzrtree_stub.size() << endl;
0800     }else{
0801       break;
0802     }
0803     xztrkmap.clear();
0804   }
0805   if( yztrkmap.size()>0){
0806     int size_before = yzrtree_stub.size();
0807     //if(Verbosity()>0)
0808     cout << "mapend: " << yztrkmap.rbegin()->first << " | " << (yztrkmap.rbegin()->second).first << " | "  << yztrkmap.rbegin()->second.second << endl;
0809     // for (const auto& [key, value] : trkmap)
0810     //  std::cout <<" ev: " << _nevent << '[' << key << "] = " << value.first << " | " << value.second << std::endl;
0811     //remove stubs in box
0812     float rmint = yztrkmap.rbegin()->second.first;
0813     float rmsl = yztrkmap.rbegin()->second.second;
0814     vector<pointKey> rmcand;
0815     yzrtree_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));
0816     for(vector<pointKey>::iterator rmstub = rmcand.begin();rmstub!=rmcand.end();++rmstub){
0817       float rmpint = rmstub->first.get<0>();
0818       float rmpsl = rmstub->first.get<1>();
0819       if(Verbosity()>0) cout<< "    rm " <<  " int: " << rmpint 
0820                 << " sl: " << rmpsl 
0821                 << endl;
0822       
0823       //auto rmpoint = std::make_pair(point(rmpint , rmpsl, 0.0), 0);
0824       //rtree_stub.remove(rmpoint);
0825       //    pointKey match;
0826       yzrtree_stub.remove(*rmstub);
0827       /*    if (!rtree_stub.query(bgi::nearest(rmpoint, 1), &match))
0828         continue;
0829       */
0830       /*
0831     if (bg::distance(rmpoint, match.first) <= 0.5) {
0832     rtree_stub.remove(match);
0833     
0834       */
0835     }
0836     int size_after = yzrtree_stub.size();
0837     if(size_before == size_after) {break;}
0838     yzouttrkmap[nout_yz++] = std::make_pair(rmint,rmsl);
0839     cout << " tree sixe after remove: " << yzrtree_stub.size() << endl;
0840   }else{
0841     break;
0842   }
0843   yztrkmap.clear();
0844   }
0845   
0846   for (const auto& [key, value] : xzouttrkmap){
0847     //if(Verbosity()>0) 
0848     std::cout <<" xz cand out ev: " << _nevent << '[' << key << "] = " << value.first << " | " << value.second << std::endl;
0849   }
0850 
0851   for (const auto& [key, value] : yzouttrkmap){
0852     //if(Verbosity()>0) 
0853     std::cout <<" yz cand out ev: " << _nevent << '[' << key << "] = " << value.first << " | " << value.second << std::endl;
0854   }
0855   
0856   for (const auto& [xzkey, xzvalue] : xzouttrkmap){
0857     std::cout <<" xz cand out ev: " << _nevent << '[' << xzkey << "] = " << xzvalue.first << " | " << xzvalue.second << std::endl;
0858 
0859     //Loop over clusters and pick up clusters compatible with line parameters
0860 
0861     float xzint = -1*xzvalue.first;
0862     float xzsl = -1*xzvalue.second;
0863     
0864     for (const auto& [yzkey, yzvalue] : yzouttrkmap){
0865       std::cout <<" yyz cand out ev: " << _nevent << '[' << xzkey << "] = " << yzvalue.first << " | " << xzvalue.second << std::endl;
0866       
0867       //Loop over clusters and pick up clusters compatible with line parameters
0868       
0869       float yzint = -1*yzvalue.first;
0870       float yzsl = -1*yzvalue.second;
0871       
0872       // The shortest distance of a point from a circle is along the radial; line from the circle center to the point
0873       /*
0874       double xza = -xzsl;//-slope;
0875       double xzb = -1.0;
0876       double xzc = -tint;//-intercept;
0877       double xzr = 100;
0878       double xzx0 = -xza*xzc/(xza*xza+xzb*xzb), xzy0 = -xzb*xzc/(xza*xza+xzb*xzb);
0879       double xzax = NAN, xzay = NAN, xzbx = NAN, xzby = NAN;
0880       if ((xzc*xzc) > (xzr*xzr*(xza*xzxa+xzb*xzb)+0.00000001)){
0881     if(Verbosity()>0) cout << "no points " << endl;
0882       }
0883       else if (abs (xzc*xzc - xzr*xzr*(xza*xza+xzb*xzxb)) < 0.0000001) {
0884       if(Verbosity()>0){ puts ("1 point");
0885     cout << xzx0 << ' ' << xzy0 << endl;
0886       }
0887       }
0888       else {
0889     double xzd = xzr*xzr - xzc*xzc/(xza*xza+xzb*xzb);
0890     double xzmult = sqrt (xzd / (xza*xza+xzb*xzb));
0891     xzax = xzx0 + xzb * xzmult;
0892     xzbx = xzx0 - xzb * xzmult;
0893     xzay = xzy0 - xza * xzmult;
0894     xzby = xzy0 + xza * xzmult;
0895     if(Verbosity()>0){puts ("2 points");
0896       cout << ax << ' ' << ay << '\n' << bx << ' ' << by << endl;
0897     }
0898       }
0899 
0900     float dist = 4;
0901     float alpha = atan(tsl);
0902     float offax = -1*dist*sin(alpha);
0903     float offay = -1*dist*cos(alpha);
0904     float offbx = 1*dist*sin(alpha);
0905     float offby = 1*dist*cos(alpha);
0906     vector<pointKey> trkclusters;
0907       */
0908       //    float trkdist = 2;
0909       /* float zstart = -110;
0910     float zend = 110;
0911     float xstart = xzint + zstart*xzsl; 
0912     float xend = xzint + zend*xzsl; 
0913     float ystart = yzint + zstart*yzsl; 
0914     float yend = yzint + zend*yzsl; */
0915     /*
0916     if(Verbosity()>0){ cout << ax << ' ' << ay << ' ' << bx << ' ' << by << endl;
0917       cout << "b1 = new TLine(" << ax +offax << ',' << ay +offay << ',' << bx -offbx << ',' << by - offby<<")" << endl;
0918       cout << "b2 = new TLine(" << ax -offax << ',' << ay -offay << ',' << bx +offbx << ',' << by + offby<<")" << endl;
0919       cout << "b1->Draw(\"same\")" << endl;
0920       cout << "b2->Draw(\"same\")" << endl;
0921       }*/
0922     vector<pointKey> trkclusters;
0923     vector<pointKey> lineclusters;
0924     rtree.query(bgi::intersects(box(point(-80,-80,-110),point(80,80,110))),std::back_inserter(lineclusters));
0925     // if(Verbosity()>0) 
0926     cout << "number line clus is " << lineclusters.size() << endl;
0927     
0928     for(vector<pointKey>::iterator clustertrk = lineclusters.begin();clustertrk!=lineclusters.end();++clustertrk){
0929       
0930       float ptx = clustertrk->first.get<0>();
0931       float pty = clustertrk->first.get<1>();
0932       float ptz = clustertrk->first.get<2>();
0933       
0934       float refx = xzint + ptz*xzsl;
0935       float refy = yzint + ptz*yzsl;
0936       float resx =  abs(ptx-refx);
0937       float resy =  abs(pty-refy);
0938       //      if(Verbosity()>0)
0939       if((resx<2)&&(resy<2)){
0940     cout << " x: " << ptx << " | " << refx << " y: " << pty << " | " << refy << " resx " << resx << " resy:" << resy << endl;
0941     trkclusters.push_back(*clustertrk);
0942       }
0943     }
0944     int nclus = 0;
0945     int nclus0 = 0;
0946     int nclus1 = 0;
0947     std::vector<float> dedxlist;
0948 
0949     cout << "number trk clus is " << trkclusters.size() << endl;
0950     if(trkclusters.size()>=_min_nclusters){
0951       cout << "setting keep true: " << trkclusters.size() << " > " << _min_nclusters << endl;
0952       for (const auto& trkclu : trkclusters){
0953     //      std::cout <<" yyz cand out ev: " << _nevent << '[' << xzkey << "] = " << yzvalue.first << " | " << xzvalue.second << std::endl;
0954     // for(vector<pointKey>::iterator clustertrk = trkclusters.begin();clustertrk!=trkclusters.end();++clustertrk){
0955     float tptx = trkclu.first.get<0>();
0956     float tpty = trkclu.first.get<1>();
0957     float tptz = trkclu.first.get<2>();
0958     //      float ptx = clustertrk->first.get<0>();
0959     //float pty = clustertrk->first.get<1>();
0960     // float ptz = clustertrk->first.get<2>();
0961      nclus++;
0962      if(tptz<0){
0963        nclus0 = 0;
0964      }
0965      if(tptz>0){
0966        nclus1 = 0;
0967      }
0968      float adc = trkclu.second;
0969      
0970      float thick = 1.098;
0971      float r  = sqrt(tptx*tptx+tpty*tpty);
0972      if(r>30&&r<=40.5){
0973        thick = 0.57;
0974      }else if(r>40.5&&r<60){
0975        thick = 1.02;
0976      } 
0977      dedxlist.push_back(adc/thick);
0978      //calc dedx
0979      //_ntp_clutrk->Fill(_nevent,tptx , tpty, tptz);
0980       }
0981       sort(dedxlist.begin(), dedxlist.end());
0982       int trunc_min = 0;
0983       int trunc_max = (int)dedxlist.size()*0.7;
0984       float sumdedx = 0;
0985       int ndedx = 0;
0986       for(int j = trunc_min; j<=trunc_max;j++){
0987     sumdedx+=dedxlist.at(j);
0988     ndedx++;
0989       }
0990       sumdedx/=ndedx;
0991       TVector2 vec(xzint,yzint);
0992       // ev:nclus:nclus0:nclus1:r:phi:dedx:intxz:intyz:intrz:slxz:slyz:slrz
0993       float r = sqrt(xzint*xzint+ yzint*yzint);
0994       float xp = xzint+xzsl*100;
0995       float xn = xzint+xzsl*-100;
0996       float yp = yzint+yzsl*100;
0997       float yn = yzint+yzsl*-100;
0998       float rp = sqrt(xp*xp+yp*yp);
0999       float rn = sqrt(xn*xn+yn*yn);
1000       float phi = vec.Phi();
1001       float dr = rp-rn;
1002       
1003       _ntp_trk->Fill(_nevent,ntrack,nclus,nclus0,nclus1,r,phi,sumdedx,xzint,yzint,r,xzsl,yzsl,dr/200);
1004 
1005       for (const auto& trkclu : trkclusters){
1006     float tptx = trkclu.first.get<0>();
1007     float tpty = trkclu.first.get<1>();
1008     float tptz = trkclu.first.get<2>();
1009     float ftX[20] = {0};
1010     int nft = 0;
1011     ftX[nft++] = _nevent;
1012     ftX[nft++] =ntrack;
1013     ftX[nft++] =tptx;
1014     ftX[nft++] =tpty;
1015     ftX[nft++] = tptz;
1016     ftX[nft++] =nclus;
1017     ftX[nft++] =nclus0;
1018     ftX[nft++] =nclus1;
1019     ftX[nft++] =r;
1020     ftX[nft++] =phi;
1021     ftX[nft++] =sumdedx;
1022     ftX[nft++] =xzint;
1023     ftX[nft++] =yzint;
1024     ftX[nft++] =r;
1025     ftX[nft++] =xzsl;
1026     ftX[nft++] =yzsl;
1027     ftX[nft++] =dr/200;
1028     
1029     _ntp_clutrk->Fill(ftX);
1030       }
1031       ntrack++;
1032       //      keep_event = true;
1033     }
1034     /*
1035     //Assemble tracks
1036     if(_create_tracks){
1037       if(trkclusters.size()>=20){
1038     auto trackseed = std::make_unique<TrackSeed_v2>();
1039     for(const auto& cluskeys:trkclusters){
1040       //      for(vector<pointKey>::iterator cluskeys =trkclusters.begin();cluskeys!=trkclusters.end();trkclusters++){
1041       
1042       trackseed->insert_cluster_key(cluskeys.second);
1043     }
1044     //    clean_chains.push_back(trackseed);
1045     // auto pseed = std::make_unique<TrackSeed_v2>(trackseed);
1046     m_seedContainer->insert(trackseed.get()); 
1047     cout << "number trk keys is " << trackseed->size_cluster_keys() << endl;
1048     numberofseeds++;
1049       }
1050     */
1051     }
1052   }
1053 
1054   
1055   cout << "number of seeds is " << ntrack << endl;
1056 
1057   /*  if(!keep_event){
1058     cout << " ABORT !keep " << endl;
1059     return Fun4AllReturnCodes::ABORTEVENT;
1060   }
1061   */
1062 #endif
1063   return Fun4AllReturnCodes::EVENT_OK;
1064 }
1065 
1066 
1067 int PHStreakFinder::Setup(PHCompositeNode* topNode)
1068 {
1069   cout << "Called Setup" << endl;
1070   cout << "topNode:" << topNode << endl;
1071   // PHTrackSeeding::Setup(topNode);
1072   //  int ret = GetNodes(topNode);
1073   //return ret;
1074   GetNodes(topNode);
1075   return Fun4AllReturnCodes::EVENT_OK;
1076 }
1077 
1078 int PHStreakFinder::End(PHCompositeNode* topNode)
1079 {
1080   if(topNode==nullptr){
1081     std::cout << PHWHERE << "No topNode, bailing."
1082           << std::endl;
1083     return Fun4AllReturnCodes::ABORTEVENT;
1084   }
1085     
1086     _tfile->cd();
1087     _ntp_clu->Write();
1088     _ntp_stub->Write();
1089     _ntp_trk->Write();
1090     _ntp_clutrk->Write();
1091     _tfile->Close();
1092     cout << "Called End " << endl;
1093  
1094   return Fun4AllReturnCodes::EVENT_OK;
1095 }
1096 
1097
[ 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!