packages/tpc/TpcClusterCleaner.cc

0001 #include "TpcClusterCleaner.h"
0002
0003 /// Tracking includes
0004
0005 #include <trackbase/TpcDefs.h>
0006 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/TrkrDefs.h>  // for cluskey, getLayer, TrkrId
0009 #include <trackbase/TrkrHitSet.h>
0010 #include <trackbase/TrkrHitSetContainer.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>
0016
0017 #include <TMatrixFfwd.h>    // for TMatrixF
0018 #include <TMatrixT.h>       // for TMatrixT, ope...
0019 #include <TMatrixTUtils.h>  // for TMatrixTRow
0020
0021 #include <cmath>     // for sqrt, fabs, atan2, cos
0022 #include <iostream>  // for operator<<, basic_ostream
0023 #include <map>       // for map
0024 #include <set>       // for _Rb_tree_const_iterator
0025 #include <utility>   // for pair, make_pair
0026
0027 //____________________________________________________________________________..
0028 TpcClusterCleaner::TpcClusterCleaner(const std::string &name)
0029   : SubsysReco(name)
0030 {
0031 }
0032
0033 //____________________________________________________________________________..
0034 TpcClusterCleaner::~TpcClusterCleaner() = default;
0035
0036 //____________________________________________________________________________..
0037 int TpcClusterCleaner::InitRun(PHCompositeNode *topNode)
0038 {
0039   int ret = GetNodes(topNode);
0040   if (ret != Fun4AllReturnCodes::EVENT_OK)
0041   {
0042     return ret;
0043   }
0044
0045   return ret;
0046 }
0047
0048 //____________________________________________________________________________..
0049 int TpcClusterCleaner::process_event(PHCompositeNode *topNode)
0050 {
0051   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0052   if (!_cluster_map)
0053   {
0054     std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0055     return Fun4AllReturnCodes::ABORTEVENT;
0056   }
0057
0058   std::set<TrkrDefs::cluskey> discard_set;
0059
0060   unsigned int count_discards = 0;
0061
0062   // loop over all TPC clusters
0063   if (Verbosity() > 0)
0064   {
0065     std::cout << std::endl
0066               << "original size of cluster map: " << _cluster_map->size() << std::endl;
0067   }
0068   for (const auto &hitsetkey : _cluster_map->getHitSetKeys(TrkrDefs::TrkrId::tpcId))
0069   {
0070     TrkrClusterContainer::ConstRange clusRange = _cluster_map->getClusters(hitsetkey);
0071
0072     for (auto clusiter = clusRange.first;
0073          clusiter != clusRange.second; ++clusiter)
0074     {
0075       TrkrDefs::cluskey cluskey = clusiter->first;
0076       TrkrCluster *cluster = clusiter->second;
0077
0078       unsigned int trkrId = TrkrDefs::getTrkrId(cluskey);
0079       unsigned int layer = TrkrDefs::getLayer(cluskey);
0080
0081       if (trkrId != TrkrDefs::tpcId)
0082       {
0083         continue;  // we want only TPC clusters
0084       }
0085
0086       if (Verbosity() > 1)
0087       {
0088         std::cout << " layer " << layer << " cluster : " << cluskey
0089                   << " ADC " << cluster->getAdc() << "       errors: r-phi " << cluster->getRPhiError() << " Z " << cluster->getZError()
0090                   << std::endl;
0091       }
0092
0093       bool discard_cluster = false;
0094
0095       // We have a TPC cluster, look for reasons to discard it
0096
0097       // errors too small
0098       // associated with very large ADC values
0099       if (cluster->getRPhiError() < _rphi_error_low_cut)
0100       {
0101         discard_cluster = true;
0102       }
0103
0104       // errors too large
0105       // associated with very small ADC values
0106       // if(cluster->getRPhiError() > _rphi_error_high_cut)
0107       // discard_cluster = true;
0108
0109       if (discard_cluster)
0110       {
0111         count_discards++;
0112         // mark it for modification
0113         discard_set.insert(cluskey);
0114         if (Verbosity() > 0)
0115         {
0116           std::cout << "                       discard cluster " << cluskey << " with ephi " << cluster->getRPhiError() << " adc " << cluster->getAdc()
0117                     << std::endl;
0118         }
0119       }
0120     }
0121   }
0122
0123   for (unsigned long iter : discard_set)
0124   {
0125     // remove bad clusters from the node tree map
0126     _cluster_map->removeCluster(iter);
0127   }
0128
0129   if (Verbosity() > 0)
0130   {
0131     std::cout << "Clusters discarded this event: " << count_discards << std::endl;
0132     std::cout << "Clusters remaining: " << _cluster_map->size() << std::endl;
0133   }
0134
0135   /*
0136   // check the map on the node tree
0137   std::cout << " Readback modified cluster map - size is: " << _cluster_map->size() << std::endl;
0138   clusRange = _cluster_map->getClusters();
0139
0140   for (auto clusiter = clusRange.first;
0141        clusiter != clusRange.second; ++clusiter)
0142     {
0143       TrkrDefs::cluskey cluskey = clusiter->first;
0144       TrkrCluster *cluster = clusiter->second;
0145
0146       unsigned int trkrId = TrkrDefs::getTrkrId(cluskey);
0147       unsigned int layer = TrkrDefs::getLayer(cluskey);
0148
0149       if(trkrId != TrkrDefs::tpcId) continue;  // we want only TPC clusters
0150
0151       if (Verbosity() >= 1)
0152         {
0153           std::cout << " cluster : " << cluskey << " layer " << layer
0154                     << " position x,y,z " << cluster->getX() << "  " << cluster->getY() << "  " << cluster->getZ()
0155                       << " ADC " << cluster->getAdc()
0156                     << std::endl;
0157           //std::cout << "       errors: r-phi " << cluster->getRPhiError() << " Z " << cluster->getZError()
0158           //          << " phi size " << cluster->getPhiSize() << " Z size " << cluster->getZSize()
0159           //          << std::endl;
0160         }
0161     }
0162   */
0163
0164   return Fun4AllReturnCodes::EVENT_OK;
0165 }
0166
0167 int TpcClusterCleaner::End(PHCompositeNode * /*topNode*/)
0168 {
0169   return Fun4AllReturnCodes::EVENT_OK;
0170 }
0171
0172 int TpcClusterCleaner::GetNodes(PHCompositeNode * /*topNode*/)
0173 {
0174   return Fun4AllReturnCodes::EVENT_OK;
0175 }
0176
0177 void TpcClusterCleaner::rotate_error(double erphi, double ez, double clusphi, double error[][3])
0178 {
0179   TMatrixF ROT(3, 3);
0180   TMatrixF ERR(3, 3);
0181   for (int i = 0; i < 3; ++i)
0182   {
0183     for (int j = 0; j < 3; ++j)
0184     {
0185       ROT[i][j] = 0;
0186       ERR[i][j] = 0;
0187     }
0188   }
0189
0190   ROT[0][0] = cos(clusphi);
0191   ROT[0][1] = -sin(clusphi);
0192   ROT[1][0] = sin(clusphi);
0193   ROT[1][1] = cos(clusphi);
0194   ROT[2][2] = 1.0;
0195
0196   ERR[1][1] = erphi * erphi;
0197   ERR[2][2] = ez * ez;
0198
0199   TMatrixF ROT_T(3, 3);
0200   ROT_T.Transpose(ROT);
0201
0202   TMatrixF COVAR_ERR(3, 3);
0203   COVAR_ERR = ROT * ERR * ROT_T;
0204
0205   for (int i = 0; i < 3; ++i)
0206   {
0207     for (int j = 0; j < 3; ++j)
0208     {
0209       error[i][j] = COVAR_ERR[i][j];
0210     }
0211   }
0212 }