packages/tpc/Tpc3DClusterizer.cc

0001 #include "Tpc3DClusterizer.h"
0002
0003 #include <trackbase/LaserCluster.h>
0004 #include <trackbase/LaserClusterContainer.h>
0005 #include <trackbase/LaserClusterContainerv1.h>
0006 #include <trackbase/LaserClusterv1.h>
0007 #include <trackbase/TpcDefs.h>
0008 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0009 #include <trackbase/TrkrHit.h>
0010 #include <trackbase/TrkrHitSet.h>
0011 #include <trackbase/TrkrHitSetContainer.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>  // for SubsysReco
0015
0016 #include <g4detectors/PHG4TpcCylinderGeom.h>
0017 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>  // for PHIODataNode
0021 #include <phool/PHNode.h>        // for PHNode
0022 #include <phool/PHNodeIterator.h>
0023 #include <phool/PHObject.h>  // for PHObject
0024 #include <phool/PHTimer.h>
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>  // for PHWHERE
0027 #include <phool/recoConsts.h>
0028
0029 #include <Acts/Definitions/Units.hpp>
0030 #include <Acts/Surfaces/Surface.hpp>
0031
0032 #include <TF1.h>
0033 #include <TFile.h>
0034
0035 #include <boost/geometry.hpp>
0036 #include <boost/geometry/geometries/box.hpp>
0037 #include <boost/geometry/geometries/point.hpp>
0038 #include <boost/geometry/index/rtree.hpp>
0039
0040 #include <array>
0041 #include <cmath>  // for sqrt, cos, sin
0042 #include <iostream>
0043 #include <limits>
0044 #include <map>  // for _Rb_tree_cons...
0045 #include <string>
0046 #include <utility>  // for pair
0047 #include <vector>
0048
0049 namespace bg = boost::geometry;
0050 namespace bgi = boost::geometry::index;
0051
0052 using point = bg::model::point<float, 3, bg::cs::cartesian>;
0053 using box = bg::model::box<point>;
0054 using specHitKey = std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey>;
0055 using pointKeyLaser = std::pair<point, specHitKey>;
0056
0057 Tpc3DClusterizer::Tpc3DClusterizer(const std::string &name)
0058   : SubsysReco(name)
0059 {
0060 }
0061
0062 int Tpc3DClusterizer::InitRun(PHCompositeNode *topNode)
0063 {
0064   PHNodeIterator iter(topNode);
0065
0066   // Looking for the DST node
0067   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0068   if (!dstNode)
0069   {
0070     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0071     return Fun4AllReturnCodes::ABORTRUN;
0072   }
0073
0074   // Create the Cluster node if required
0075   auto laserclusters = findNode::getClass<LaserClusterContainer>(dstNode, "LASER_CLUSTER");
0076   if (!laserclusters)
0077   {
0078     PHNodeIterator dstiter(dstNode);
0079     PHCompositeNode *DetNode =
0080         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0081     if (!DetNode)
0082     {
0083       DetNode = new PHCompositeNode("TRKR");
0084       dstNode->addNode(DetNode);
0085     }
0086
0087     laserclusters = new LaserClusterContainerv1;
0088     PHIODataNode<PHObject> *LaserClusterContainerNode =
0089         new PHIODataNode<PHObject>(laserclusters, "LASER_CLUSTER", "PHObject");
0090     DetNode->addNode(LaserClusterContainerNode);
0091   }
0092
0093   if (m_debug)
0094   {
0095     m_debugFile = new TFile(m_debugFileName.c_str(), "RECREATE");
0096   }
0097
0098   m_itHist_0 = new TH1I("m_itHist_0", "side 0;it", 360, -0.5, 359.5);
0099   m_itHist_1 = new TH1I("m_itHist_1", "side 1;it", 360, -0.5, 359.5);
0100
0101   if (m_debug)
0102   {
0103     m_clusterTree = new TTree("clusterTree", "clusterTree");
0104     m_clusterTree->Branch("event", &m_event);
0105     m_clusterTree->Branch("clusters", &m_eventClusters);
0106     m_clusterTree->Branch("itHist_0", &m_itHist_0);
0107     m_clusterTree->Branch("itHist_1", &m_itHist_1);
0108     m_clusterTree->Branch("nClusters", &m_nClus);
0109     m_clusterTree->Branch("time_search", &time_search);
0110     m_clusterTree->Branch("time_clus", &time_clus);
0111     m_clusterTree->Branch("time_erase", &time_erase);
0112     m_clusterTree->Branch("time_all", &time_all);
0113   }
0114
0115   if (m_output){
0116     m_outputFile = new TFile(m_outputFileName.c_str(), "RECREATE");
0117
0118     m_clusterNT = new TNtuple("clus3D", "clus3D","event:seed:x:y:z:r:phi:phibin:tbin:adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize");
0119   }
0120
0121   m_geom_container =
0122       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0123   if (!m_geom_container)
0124   {
0125     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0126     return Fun4AllReturnCodes::ABORTRUN;
0127   }
0128
0129   AdcClockPeriod = m_geom_container->GetFirstLayerCellGeom()->get_zstep();
0130
0131   m_tdriftmax = AdcClockPeriod * NZBinsSide;
0132
0133   t_all = std::make_unique<PHTimer>("t_all");
0134   t_all->stop();
0135   t_search = std::make_unique<PHTimer>("t_search");
0136   t_search->stop();
0137   t_clus = std::make_unique<PHTimer>("t_clus");
0138   t_clus->stop();
0139   t_erase = std::make_unique<PHTimer>("t_erase");
0140   t_erase->stop();
0141
0142   return Fun4AllReturnCodes::EVENT_OK;
0143 }
0144
0145 int Tpc3DClusterizer::process_event(PHCompositeNode *topNode)
0146 {
0147   ++m_event;
0148   recoConsts* rc = recoConsts::instance();
0149   if (rc->FlagExist("RANDOMSEED")){
0150     m_seed = (int)rc->get_IntFlag("RANDOMSEED");
0151   }  else {
0152     m_seed = std::numeric_limits<int>::quiet_NaN();
0153   }
0154
0155   std::cout << "Tpc3DClusterizer::process_event working on event " << m_event << std::endl;
0156
0157   PHNodeIterator iter(topNode);
0158   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0159   if (!dstNode)
0160   {
0161     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0162     return Fun4AllReturnCodes::ABORTRUN;
0163   }
0164   // get node containing the digitized hits
0165   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0166   if (!m_hits){
0167     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0168     return Fun4AllReturnCodes::ABORTRUN;
0169   }
0170
0171   // get node for clusters
0172   m_clusterlist = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0173   if (!m_clusterlist)
0174   {
0175     std::cout << PHWHERE << " ERROR: Can't find LASER_CLUSTER." << std::endl;
0176     return Fun4AllReturnCodes::ABORTRUN;
0177   }
0178
0179   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
0180                                                  "ActsGeometry");
0181   if (!m_tGeometry)
0182   {
0183     std::cout << PHWHERE
0184               << "ActsGeometry not found on node tree. Exiting"
0185               << std::endl;
0186     return Fun4AllReturnCodes::ABORTRUN;
0187   }
0188
0189   TrkrHitSetContainer::ConstRange hitsetrange;
0190   hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0191
0192   bgi::rtree<pointKeyLaser, bgi::quadratic<16>> rtree;
0193   bgi::rtree<pointKeyLaser, bgi::quadratic<16>> rtree_reject;
0194   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0195        hitsetitr != hitsetrange.second;
0196        ++hitsetitr){
0197     TrkrHitSet *hitset = hitsetitr->second;
0198     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0199     int side = TpcDefs::getSide(hitsetitr->first);
0200     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0201     TrkrDefs::hitsetkey hitsetKey = TpcDefs::genHitSetKey(layer, sector, side);
0202
0203     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0204     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0205      hitr != hitrangei.second;
0206      ++hitr){
0207       int iphi = TpcDefs::getPad(hitr->first);
0208       int it = TpcDefs::getTBin(hitr->first);
0209       //      std::cout << " iphi: " << iphi << " it: " << it << std::endl;
0210       float_t fadc = (hitr->second->getAdc());// - m_pedestal;  // proper int rounding +0.5
0211       unsigned short adc = 0;
0212       if (fadc > 0){
0213     adc = (unsigned short) fadc;
0214       }
0215       if (adc <= 0){
0216     continue;
0217       }
0218
0219       std::vector<pointKeyLaser> testduplicate;
0220       rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001),
0221                       point(layer + 0.001, iphi + 0.001, it + 0.001))),
0222           std::back_inserter(testduplicate));
0223       if (!testduplicate.empty()){
0224     testduplicate.clear();
0225     continue;
0226       }
0227
0228       TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it);
0229
0230       auto spechitkey = std::make_pair(hitKey, hitsetKey);
0231       rtree_reject.insert(std::make_pair(point(1.0 * layer, 1.0 * iphi, 1.0 * it), spechitkey));
0232     }
0233   }
0234
0235   std::multimap<unsigned int, std::pair<std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey>, std::array<int, 3>>> adcMap;
0236   //  std::cout << "n hitsets: " << std::distance(hitsetrange.first,hitsetrange.second)
0237   //          << std::endl;
0238   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0239        hitsetitr != hitsetrange.second;
0240        ++hitsetitr){
0241     TrkrHitSet *hitset = hitsetitr->second;
0242     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0243     int side = TpcDefs::getSide(hitsetitr->first);
0244     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0245     //PHG4TpcCylinderGeom *layergeom = m_geom_container->GetLayerCellGeom(layer);
0246     // double r = layergeom->get_radius();
0247
0248     TrkrDefs::hitsetkey hitsetKey = TpcDefs::genHitSetKey(layer, sector, side);
0249
0250     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0251     // std::cout << "n hits in set: " << hitset->size()
0252     //         << std::endl;
0253     // int nhits = 0;
0254     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0255      hitr != hitrangei.second;
0256      ++hitr){
0257       int iphi = TpcDefs::getPad(hitr->first);
0258       int it = TpcDefs::getTBin(hitr->first);
0259       // std::cout << " iphi: " << iphi << " it: " << it << std::endl;
0260       float_t fadc = (hitr->second->getAdc());// - m_pedestal;  // proper int rounding +0.5
0261       unsigned short adc = 0;
0262       // std::cout << " nhit: " << nhits++ << "adc: " << fadc << " phi: " << iphi << " it: " << it << std::endl;
0263       if (fadc > 0){
0264     adc = (unsigned short) fadc;
0265       }
0266       if (adc <= 0){
0267     continue;
0268       }
0269       if(layer>=7+32){
0270     //if(side==1)continue;
0271     if(abs(iphi-0)<=2) continue;
0272     if(abs(iphi-191)<=2) continue;
0273     if(abs(iphi-206)<=1) continue;
0274     if(abs(iphi-383)<=2) continue;
0275     if(abs(iphi-576)<=2) continue;
0276     if(abs(iphi-767)<=2) continue;
0277     if(abs(iphi-960)<=2) continue;
0278     if(abs(iphi-1522)<=2) continue;
0279     if(abs(iphi-1344)<=2) continue;
0280     if(abs(iphi-1536)<=2) continue;
0281     if(abs(iphi-1728)<=2) continue;
0282     if(abs(iphi-1920)<=2) continue;
0283     if(abs(iphi-2111)<=2) continue;
0284     if(abs(iphi-2303)<=2) continue;
0285       }
0286
0287       /*
0288       double phi = layergeom->get_phi(iphi);
0289       double m_sampa_tbias = 39.6;
0290       double zdriftlength = (layergeom->get_zcenter(it)+ m_sampa_tbias) * m_tGeometry->get_drift_velocity();
0291
0292       float x = r * cos(phi);
0293       float y = r * sin(phi);
0294       float z = m_tdriftmax * m_tGeometry->get_drift_velocity() - zdriftlength;
0295       if (side == 0){
0296     z = -z;
0297     it = -it;
0298       }
0299       */
0300       std::array<int, 3> coords = {(int) layer, iphi, it};
0301
0302       std::vector<pointKeyLaser> testduplicate;
0303       rtree.query(bgi::intersects(box(point(layer - 0.001, iphi - 0.001, it - 0.001),
0304                       point(layer + 0.001, iphi + 0.001, it + 0.001))),
0305           std::back_inserter(testduplicate));
0306       if (!testduplicate.empty()){
0307     testduplicate.clear();
0308     continue;
0309       }
0310
0311       //test for isolated hit
0312       std::vector<pointKeyLaser> testisolated;
0313       rtree_reject.query(bgi::intersects(box(point(layer - 1.001, iphi - 1.001, it - 1.001),
0314                       point(layer + 1.001, iphi + 1.001, it + 1.001))),
0315           std::back_inserter(testisolated));
0316       if(testisolated.size()==1){
0317     //testisolated.clear();
0318     continue;
0319       }
0320
0321       TrkrDefs::hitkey hitKey = TpcDefs::genHitKey(iphi, it);
0322
0323       auto spechitkey = std::make_pair(hitKey, hitsetKey);
0324       auto keyCoords = std::make_pair(spechitkey, coords);
0325       adcMap.insert(std::make_pair(adc, keyCoords));
0326       // std::cout << "inserting " << " l: " << layer << " iphi: " <<iphi << " t: " << it << std::endl;
0327       rtree.insert(std::make_pair(point(1.0 * layer, 1.0 * iphi, 1.0 * it), spechitkey));
0328     }
0329   }
0330
0331   if (Verbosity() > 1){
0332     std::cout << "finished looping over hits" << std::endl;
0333     std::cout << "map size: " << adcMap.size() << std::endl;
0334     std::cout << "rtree size: " << rtree.size() << std::endl;
0335   }
0336
0337   // done filling rTree
0338
0339   t_all->restart();
0340
0341   while (adcMap.size() > 0){
0342     auto iterKey = adcMap.rbegin();
0343     if (iterKey == adcMap.rend()){
0344       break;
0345     }
0346
0347     auto coords = iterKey->second.second;
0348
0349     int layer = coords[0];
0350     int iphi = coords[1];
0351     int it = coords[2];
0352
0353     int layerMax = layer + 1;
0354     if (layer == 22 || layer == 38 || layer == 54){
0355       layerMax = layer;
0356     }
0357     int layerMin = layer - 1;
0358     if (layer == 7 || layer == 23 || layer == 39){
0359       layerMin = layer;
0360     }
0361
0362     std::vector<pointKeyLaser> clusHits;
0363
0364     t_search->restart();
0365     rtree.query(bgi::intersects(box(point(layerMin, iphi - 2, it - 5), point(layerMax, iphi + 2, it + 5))), std::back_inserter(clusHits));
0366     t_search->stop();
0367
0368     t_clus->restart();
0369     calc_cluster_parameter(clusHits, adcMap);
0370     t_clus->stop();
0371
0372     t_erase->restart();
0373     remove_hits(clusHits, rtree, adcMap);
0374     t_erase->stop();
0375
0376     clusHits.clear();
0377   }
0378
0379   if (m_debug){
0380     m_nClus = (int) m_eventClusters.size();
0381   }
0382   t_all->stop();
0383
0384   if (m_debug){
0385     time_search = t_search->get_accumulated_time() / 1000.;
0386     time_clus = t_clus->get_accumulated_time() / 1000.;
0387     time_erase = t_erase->get_accumulated_time() / 1000.;
0388     time_all = t_all->get_accumulated_time() / 1000.;
0389
0390     m_clusterTree->Fill();
0391   }
0392
0393   if (Verbosity()){
0394     std::cout << "rtree search time: " << t_search->get_accumulated_time() / 1000. << " sec" << std::endl;
0395     std::cout << "clustering time: " << t_clus->get_accumulated_time() / 1000. << " sec" << std::endl;
0396     std::cout << "erasing time: " << t_erase->get_accumulated_time() / 1000. << " sec" << std::endl;
0397     std::cout << "total time: " << t_all->get_accumulated_time() / 1000. << " sec" << std::endl;
0398   }
0399
0400   //  std::cout << " not enough clusters in _nevent: " << m_clusterlist->size() << std::endl;
0401   return Fun4AllReturnCodes::EVENT_OK;
0402 }
0403
0404 int Tpc3DClusterizer::ResetEvent(PHCompositeNode * /*topNode*/){
0405   m_itHist_0->Reset();
0406   m_itHist_1->Reset();
0407
0408   if (m_debug)
0409     {
0410       m_tHist_0->Reset();
0411       m_tHist_1->Reset();
0412
0413       m_eventClusters.clear();
0414     }
0415
0416   return Fun4AllReturnCodes::EVENT_OK;
0417 }
0418
0419 int Tpc3DClusterizer::End(PHCompositeNode * /*topNode*/)
0420 {
0421   if (m_debug){
0422     m_debugFile->cd();
0423     m_clusterTree->Write();
0424     m_debugFile->Close();
0425   }
0426
0427   if (m_output){
0428     m_outputFile->cd();
0429     m_clusterNT->Write();
0430     m_outputFile->Close();
0431   }
0432   return Fun4AllReturnCodes::EVENT_OK;
0433 }
0434
0435 void Tpc3DClusterizer::calc_cluster_parameter(std::vector<pointKeyLaser> &clusHits, std::multimap<unsigned int, std::pair<std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey>, std::array<int, 3>>> &adcMap)
0436 {
0437   //std::cout << "nu clus" << std::endl;
0438   double rSum = 0.0;
0439   double phiSum = 0.0;
0440   double tSum = 0.0;
0441
0442   double layerSum = 0.0;
0443   double iphiSum = 0.0;
0444   double itSum = 0.0;
0445
0446   double adcSum = 0.0;
0447
0448   double maxAdc = 0.0;
0449   TrkrDefs::hitsetkey maxKey = 0;
0450
0451   unsigned int nHits = clusHits.size();
0452   int iphimin = 6666, iphimax = -1;
0453   int ilaymin = 6666, ilaymax = -1;
0454   float itmin = 66666666.6, itmax = -6666666666.6;
0455   auto *clus = new LaserClusterv1;
0456
0457   for (auto &clusHit : clusHits)
0458   {
0459     float coords[3] = {clusHit.first.get<0>(), clusHit.first.get<1>(), clusHit.first.get<2>()};
0460     std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey> spechitkey = clusHit.second;
0461
0462     int side = TpcDefs::getSide(spechitkey.second);
0463     // unsigned int sector= TpcDefs::getSectorId(spechitkey.second);
0464
0465     PHG4TpcCylinderGeom *layergeom = m_geom_container->GetLayerCellGeom((int) coords[0]);
0466
0467     double r = layergeom->get_radius();
0468     double phi = layergeom->get_phi(coords[1], side);
0469     double t = layergeom->get_zcenter(fabs(coords[2]));
0470     int tbin = coords[2];
0471     int lay = coords[0];//TrkrDefs::getLayer(spechitkey.second);
0472     double hitzdriftlength = t * m_tGeometry->get_drift_velocity();
0473     double hitZ = m_tdriftmax * m_tGeometry->get_drift_velocity() - hitzdriftlength;
0474     /*std::cout << " lay: " << lay
0475           << " phi: " << phi
0476           << " t: " << t
0477           << " side: " << side
0478           << std::endl;
0479     */
0480     if(phi<iphimin){iphimin = phi;}
0481     if(phi>iphimax){iphimax = phi;}
0482     if(lay<ilaymin){ilaymin = lay;}
0483     if(lay>ilaymax){ilaymax = lay;}
0484     if(tbin<itmin){itmin = tbin;}
0485     if(tbin>itmax){itmax = tbin;}
0486
0487     for (auto &iterKey : adcMap)
0488     {
0489       if (iterKey.second.first == spechitkey)
0490       {
0491         double adc = iterKey.first;
0492
0493         clus->addHit();
0494         clus->setHitLayer(clus->getNhits() - 1, coords[0]);
0495         clus->setHitIPhi(clus->getNhits() - 1, coords[1]);
0496         clus->setHitIT(clus->getNhits() - 1, coords[2]);
0497         clus->setHitX(clus->getNhits() - 1, r * cos(phi));
0498         clus->setHitY(clus->getNhits() - 1, r * sin(phi));
0499         clus->setHitZ(clus->getNhits() - 1, hitZ);
0500         clus->setHitAdc(clus->getNhits() - 1, (float) adc);
0501
0502         rSum += r * adc;
0503         phiSum += phi * adc;
0504         tSum += t * adc;
0505
0506         layerSum += coords[0] * adc;
0507         iphiSum += coords[1] * adc;
0508         itSum += coords[2] * adc;
0509
0510         adcSum += adc;
0511
0512         if (adc > maxAdc)
0513         {
0514           maxAdc = adc;
0515           maxKey = spechitkey.second;
0516         }
0517
0518         break;
0519       }
0520     }
0521   }
0522
0523   if (nHits == 0)
0524   {
0525     std::cout << "no hits"<< std::endl;
0526     return;
0527   }
0528
0529   double clusR = rSum / adcSum;
0530   double clusPhi = phiSum / adcSum;
0531   double clusiPhi = iphiSum / adcSum;
0532   double clusT = tSum / adcSum;
0533   double zdriftlength = clusT * m_tGeometry->get_drift_velocity();
0534
0535   double clusX = clusR * cos(clusPhi);
0536   double clusY = clusR * sin(clusPhi);
0537   double clusZ = m_tdriftmax * m_tGeometry->get_drift_velocity() - zdriftlength;
0538   int maxside = TpcDefs::getSide(maxKey);
0539   int maxsector = TpcDefs::getSectorId(maxKey);
0540
0541   if (maxside == 0)
0542   {
0543     clusZ = -clusZ;
0544     for (int i = 0; i < (int) clus->getNhits(); i++)
0545     {
0546       clus->setHitZ(i, -1 * clus->getHitZ(i));
0547     }
0548   }
0549
0550   clus->setAdc(adcSum);
0551   clus->setX(clusX);
0552   clus->setY(clusY);
0553   clus->setZ(clusZ);
0554   clus->setLayer(layerSum / adcSum);
0555   clus->setIPhi(iphiSum / adcSum);
0556   clus->setIT(itSum / adcSum);
0557   int phisize =  iphimax - iphimin + 1;
0558   int lsize =  ilaymax - ilaymin + 1;
0559   int tsize = itmax - itmin +1;
0560   if (m_debug)
0561   {
0562     m_currentCluster = (LaserCluster *) clus->CloneMe();
0563     m_eventClusters.push_back((LaserCluster *) m_currentCluster->CloneMe());
0564   }
0565   //  if(nHits>1&&tsize>5){
0566   if(nHits>=1){
0567     const auto ckey = TrkrDefs::genClusKey(maxKey, m_clusterlist->size());
0568     m_clusterlist->addClusterSpecifyKey(ckey, clus);
0569   } else {
0570     delete clus;
0571   }
0572
0573
0574    //event:seed:x:y:z:r:phi:phibin:tbin::adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize
0575   //"event:seed:x:y:z:r:phi:phibin:tbin::adc:maxadc:layer:phielem:zelem:size:phisize:tsize:lsize"
0576   /*
0577   std::cout << " l size: " << lsize
0578         << " phisize : " << phisize
0579         << " tsize: " << tsize
0580         << " maxside: " << maxside
0581         << std::endl;
0582   */
0583   // if (m_output){
0584     float fX[20] = {0};
0585     int n = 0;
0586     fX[n++] = m_event;
0587     fX[n++] = m_seed;
0588     fX[n++] = clusX;
0589     fX[n++] = clusY;
0590     fX[n++] = clusZ;
0591     fX[n++] = clusR;
0592     fX[n++] = clusPhi;
0593     fX[n++] = clusiPhi;
0594     fX[n++] = clusT;
0595     fX[n++] = adcSum;
0596     fX[n++] = maxAdc;
0597     fX[n++] = (layerSum/adcSum);
0598     fX[n++] = maxsector;
0599     fX[n++] = maxside;
0600     fX[n++] = nHits;
0601     fX[n++] = phisize;
0602     fX[n++] = tsize;
0603     fX[n++] = lsize;
0604     m_clusterNT->Fill(fX);
0605     // }
0606 }
0607
0608 void Tpc3DClusterizer::remove_hits(std::vector<pointKeyLaser> &clusHits, bgi::rtree<pointKeyLaser, bgi::quadratic<16>> &rtree, std::multimap<unsigned int, std::pair<std::pair<TrkrDefs::hitkey, TrkrDefs::hitsetkey>, std::array<int, 3>>> &adcMap)
0609 {
0610   for (auto &clusHit : clusHits)
0611   {
0612     auto spechitkey = clusHit.second;
0613
0614     if(rtree.size()==0){
0615       std::cout << "not good" << std::endl;
0616     }
0617     //rtree.remove(clusHit);
0618
0619     for (auto iterAdc = adcMap.begin(); iterAdc != adcMap.end();)
0620     {
0621       if (iterAdc->second.first == spechitkey)
0622       {
0623         iterAdc = adcMap.erase(iterAdc);
0624         break;
0625       }
0626       else
0627       {
0628         ++iterAdc;
0629       }
0630     }
0631   }
0632 }