EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcClusterCleaner.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcClusterCleaner.cc
1 #include "TpcClusterCleaner.h"
2 
4 
5 #include <trackbase/TrkrCluster.h> // for TrkrCluster
6 #include <trackbase/TrkrDefs.h> // for cluskey, getLayer, TrkrId
8 #include <trackbase/TrkrHitSet.h>
10 #include <tpc/TpcDefs.h>
11 
13 
14 #include <phool/getClass.h>
15 #include <phool/phool.h>
16 
17 #include <TMatrixFfwd.h> // for TMatrixF
18 #include <TMatrixT.h> // for TMatrixT, ope...
19 #include <TMatrixTUtils.h> // for TMatrixTRow
20 
21 #include <cmath> // for sqrt, fabs, atan2, cos
22 #include <iostream> // for operator<<, basic_ostream
23 #include <map> // for map
24 #include <set> // for _Rb_tree_const_iterator
25 #include <utility> // for pair, make_pair
26 
27 //____________________________________________________________________________..
29  : SubsysReco(name)
30 {
31 
32 }
33 
34 //____________________________________________________________________________..
36 {
37 
38 }
39 
40 //____________________________________________________________________________..
42 {
43  int ret = GetNodes(topNode);
44  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
45 
46  return ret;
47 }
48 
49 //____________________________________________________________________________..
51 {
52 
53  _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
54  if (!_cluster_map)
55  {
56  std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
58  }
59  _hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
60  if (!_hitsets)
61  {
62  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
64  }
65 
66  std::set<TrkrDefs::cluskey> discard_set;
67 
68  unsigned int count_discards = 0;
69 
70  // loop over all TPC clusters
71  if(Verbosity() > 0) std::cout << std::endl << "original size of cluster map: " << _cluster_map->size() << std::endl;
73  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
74  hitsetitr != hitsetrange.second;
75  ++hitsetitr){
76  TrkrClusterContainer::ConstRange clusRange = _cluster_map->getClusters(hitsetitr->first);
78 
79  for (clusiter = clusRange.first;
80  clusiter != clusRange.second; ++clusiter)
81  {
82  TrkrDefs::cluskey cluskey = clusiter->first;
83  TrkrCluster *cluster = clusiter->second;
84 
85  unsigned int trkrId = TrkrDefs::getTrkrId(cluskey);
86  unsigned int layer = TrkrDefs::getLayer(cluskey);
87 
88  if(trkrId != TrkrDefs::tpcId) continue; // we want only TPC clusters
89 
90  if (Verbosity() > 1)
91  {
92  std::cout << " layer " << layer << " cluster : " << cluskey
93  << " ADC " << cluster->getAdc() << " errors: r-phi " << cluster->getRPhiError() << " Z " << cluster->getZError()
94  << std::endl;
95  }
96 
97  bool discard_cluster = false;
98 
99  // We have a TPC cluster, look for reasons to discard it
100 
101  // errors too small
102  // associated with very large ADC values
103  if(cluster->getRPhiError() < _rphi_error_low_cut)
104  discard_cluster = true;
105 
106  // errors too large
107  // associated with very small ADC values
108  //if(cluster->getRPhiError() > _rphi_error_high_cut)
109  //discard_cluster = true;
110 
111  if(discard_cluster)
112  {
113  count_discards++;
114  // mark it for modification
115  discard_set.insert(cluskey);
116  if(Verbosity() > 0)
117  std::cout << " discard cluster " << cluskey << " with ephi " << cluster->getRPhiError() << " adc " << cluster->getAdc()
118  << std::endl;
119  }
120  }
121  }
122 
123  for(auto iter = discard_set.begin(); iter != discard_set.end(); ++iter)
124  {
125 
126  // remove bad clusters from the node tree map
127  _cluster_map->removeCluster(*iter);
128 
129  }
130 
131  if(Verbosity() > 0)
132  {
133  std::cout << "Clusters discarded this event: " << count_discards << std::endl;
134  std::cout << "Clusters remaining: " << _cluster_map->size() << std::endl;
135  }
136 
137  /*
138  // check the map on the node tree
139  std::cout << " Readback modified cluster map - size is: " << _cluster_map->size() << std::endl;
140  clusRange = _cluster_map->getClusters();
141 
142  for (auto clusiter = clusRange.first;
143  clusiter != clusRange.second; ++clusiter)
144  {
145  TrkrDefs::cluskey cluskey = clusiter->first;
146  TrkrCluster *cluster = clusiter->second;
147 
148  unsigned int trkrId = TrkrDefs::getTrkrId(cluskey);
149  unsigned int layer = TrkrDefs::getLayer(cluskey);
150 
151  if(trkrId != TrkrDefs::tpcId) continue; // we want only TPC clusters
152 
153  if (Verbosity() >= 1)
154  {
155  std::cout << " cluster : " << cluskey << " layer " << layer
156  << " position x,y,z " << cluster->getX() << " " << cluster->getY() << " " << cluster->getZ()
157  << " ADC " << cluster->getAdc()
158  << std::endl;
159  //std::cout << " errors: r-phi " << cluster->getRPhiError() << " Z " << cluster->getZError()
160  // << " phi size " << cluster->getPhiSize() << " Z size " << cluster->getZSize()
161  // << std::endl;
162  }
163  }
164  */
165 
167 }
168 
170 {
172 }
173 
175 {
176 
178 }
179 
180 void TpcClusterCleaner::rotate_error(double erphi, double ez, double clusphi, double error[][3])
181 {
182  TMatrixF ROT(3, 3);
183  TMatrixF ERR(3,3);
184  for(int i=0;i<3;++i)
185  for(int j=0;j<3;++j)
186  {
187  ROT[i][j] = 0;
188  ERR[i][j] = 0;
189  }
190 
191  ROT[0][0] = cos(clusphi);
192  ROT[0][1] = -sin(clusphi);
193  ROT[1][0] = sin(clusphi);
194  ROT[1][1] = cos(clusphi);
195  ROT[2][2] = 1.0;
196 
197  ERR[1][1] = erphi*erphi;
198  ERR[2][2] = ez*ez;
199 
200  TMatrixF ROT_T(3, 3);
201  ROT_T.Transpose(ROT);
202 
203  TMatrixF COVAR_ERR(3, 3);
204  COVAR_ERR = ROT * ERR * ROT_T;
205 
206  for(int i=0;i<3;++i)
207  for(int j=0;j<3;++j)
208  error[i][j] = COVAR_ERR[i][j];
209 
210 }