ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawClusterBuilderFwd.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawClusterBuilderFwd.cc
1 #include "RawClusterBuilderFwd.h"
2 
3 #include "PHMakeGroups.h"
4 
5 #include <calobase/RawCluster.h>
6 #include <calobase/RawClusterContainer.h>
7 #include <calobase/RawClusterDefs.h>
8 #include <calobase/RawClusterv1.h>
9 #include <calobase/RawTower.h>
10 #include <calobase/RawTowerContainer.h>
11 #include <calobase/RawTowerDefs.h>
12 #include <calobase/RawTowerGeom.h>
13 #include <calobase/RawTowerGeomContainer.h>
14 
16 #include <fun4all/SubsysReco.h>
17 
18 #include <phool/getClass.h>
19 #include <phool/PHCompositeNode.h>
20 #include <phool/PHIODataNode.h>
21 #include <phool/PHNode.h>
22 #include <phool/PHNodeIterator.h>
23 #include <phool/PHObject.h>
24 #include <phool/phool.h>
25 
26 #include <cassert>
27 #include <cmath>
28 #include <exception>
29 #include <iostream>
30 #include <map>
31 #include <stdexcept>
32 #include <utility>
33 #include <vector>
34 
35 using namespace std;
36 
37 class twrs_fwd
38 {
39  public:
40  twrs_fwd(RawTower *);
41  virtual ~twrs_fwd() {}
42  bool is_adjacent(twrs_fwd &);
43  void set_id(const int i)
44  {
45  id = i;
46  }
47  int get_id() const
48  {
49  return id;
50  }
51  int get_j_bin() const
52  {
53  return bin_j;
54  }
55  int get_k_bin() const
56  {
57  return bin_k;
58  }
59 
60  protected:
61  int bin_j;
62  int bin_k;
64 };
65 
67  : id(-1)
68 {
69  bin_j = rt->get_bineta();
70  bin_k = rt->get_binphi();
71 }
72 
74 {
75  if (bin_j - 1 <= tower.get_j_bin() && tower.get_j_bin() <= bin_j + 1)
76  {
77  if (bin_k - 1 <= tower.get_k_bin() && tower.get_k_bin() <= bin_k + 1)
78  {
79  return true;
80  }
81  }
82 
83  return false;
84 }
85 
86 bool operator<(const twrs_fwd &a, const twrs_fwd &b)
87 {
88  if (a.get_j_bin() != b.get_j_bin())
89  {
90  return a.get_j_bin() < b.get_j_bin();
91  }
92  return a.get_k_bin() < b.get_k_bin();
93 }
94 
96  : SubsysReco(name)
97  , _clusters(nullptr)
98  , _min_tower_e(0.0)
99  , chkenergyconservation(0)
100  , detector("NONE")
101 {
102 }
103 
105 {
106  try
107  {
108  CreateNodes(topNode);
109  }
110  catch (std::exception &e)
111  {
112  std::cout << PHWHERE << ": " << e.what() << std::endl;
113  throw;
114  }
115 
117 }
118 
120 {
121  string towernodename = "TOWER_CALIB_" + detector;
122  // Grab the towers
123  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, towernodename);
124  if (!towers)
125  {
126  std::cout << PHWHERE << ": Could not find node " << towernodename << std::endl;
128  }
129  string towergeomnodename = "TOWERGEOM_" + detector;
130  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
131  if (!towergeom)
132  {
133  cout << PHWHERE << ": Could not find node " << towergeomnodename << endl;
135  }
136  // make the list of towers above threshold
137  std::vector<twrs_fwd> towerVector;
138  RawTowerContainer::ConstRange begin_end = towers->getTowers();
139  RawTowerContainer::ConstIterator itr = begin_end.first;
140  for (; itr != begin_end.second; ++itr)
141  {
142  RawTower *tower = itr->second;
143  RawTowerDefs::keytype towerid = itr->first;
144  if (tower->get_energy() > _min_tower_e)
145  {
146  twrs_fwd twr(tower);
147  twr.set_id(towerid);
148  towerVector.push_back(twr);
149  }
150  }
151 
152  // cluster the towers
153  std::multimap<int, twrs_fwd> clusteredTowers;
154  PHMakeGroups(towerVector, clusteredTowers);
155 
156  RawCluster *cluster = nullptr;
157  int last_id = -1;
158  std::multimap<int, twrs_fwd>::iterator ctitr = clusteredTowers.begin();
159  std::multimap<int, twrs_fwd>::iterator lastct = clusteredTowers.end();
160  for (; ctitr != lastct; ++ctitr)
161  {
162  int clusterid = ctitr->first;
163 
164  if (last_id != clusterid)
165  {
166  // new cluster
167  cluster = new RawClusterv1();
168  _clusters->AddCluster(cluster);
169 
170  last_id = clusterid;
171  }
172  assert(cluster);
173 
174  const twrs_fwd &tmptower = ctitr->second;
175  RawTower *rawtower = towers->getTower(tmptower.get_id());
176 
177  const double e = rawtower->get_energy();
178  cluster->addTower(rawtower->get_id(), e);
179  }
180 
181  for (const auto &cluster_pair : _clusters->getClustersMap())
182  {
183  const RawClusterDefs::keytype clusterid = cluster_pair.first;
184  RawCluster *cluster = cluster_pair.second;
185 
186  assert(cluster);
187  assert(cluster->get_id() == clusterid);
188 
189  double sum_x(0);
190  double sum_y(0);
191  double sum_z(0);
192  double sum_e(0);
193 
194  for (const auto tower_pair : cluster->get_towermap())
195  {
196  const RawTower *rawtower = towers->getTower(tower_pair.first);
197  const RawTowerGeom *rawtowergeom = towergeom->get_tower_geometry(tower_pair.first);
198 
199  assert(rawtower);
200  assert(rawtowergeom);
201  const double e = rawtower->get_energy();
202 
203  sum_e += e;
204 
205  if (e > 0)
206  {
207  sum_x += e * rawtowergeom->get_center_x();
208  sum_y += e * rawtowergeom->get_center_y();
209  sum_z += e * rawtowergeom->get_center_z();
210  }
211  } // for (const auto tower_pair : cluster->get_towermap())
212 
213  cluster->set_energy(sum_e);
214 
215  if (sum_e > 0)
216  {
217  sum_x /= sum_e;
218  sum_y /= sum_e;
219  sum_z /= sum_e;
220 
221  cluster->set_r(sqrt(sum_y * sum_y + sum_x * sum_x));
222  cluster->set_phi(atan2(sum_y, sum_x));
223  cluster->set_z(sum_z);
224  }
225 
226  if (Verbosity() > 1)
227  {
228  cout << "RawClusterBuilder constucted ";
229  cluster->identify();
230  }
231  } // for (const auto & cluster_pair : _clusters->getClustersMap())
233  {
234  double ecluster = _clusters->getTotalEdep();
235  double etower = towers->getTotalEdep();
236  if (ecluster > 0)
237  {
238  if (fabs(etower - ecluster) / ecluster > 1e-9)
239  {
240  cout << "energy conservation violation: ETower: " << etower
241  << " ECluster: " << ecluster
242  << " diff: " << etower - ecluster << endl;
243  }
244  }
245  else
246  {
247  if (etower != 0)
248  {
249  cout << "energy conservation violation: ETower: " << etower
250  << " ECluster: " << ecluster << endl;
251  }
252  }
253  }
255 }
256 
258 {
259  double sum = cluster->get_energy();
260  double phimin = 999.;
261  double phimax = -999.;
262  RawCluster::TowerConstRange begin_end = cluster->get_towers();
264  for (iter = begin_end.first; iter != begin_end.second; ++iter)
265  {
266  RawTower *tmpt = towers->getTower(iter->first);
267  RawTowerGeom *tgeo =
268  towergeom->get_tower_geometry(tmpt->get_id());
269  double phi = tgeo->get_phi();
270  if (phi > M_PI) phi = phi - 2. * M_PI; // correct the cluster phi for slat geometry which is 0-2pi (L. Xue)
271  if (phi < phimin)
272  {
273  phimin = phi;
274  }
275  if (phi > phimax)
276  {
277  phimax = phi;
278  }
279  }
280 
281  if ((phimax - phimin) < 3.) return false; // cluster is not at phi discontinuity
282 
283  float mean = 0.;
284  for (iter = begin_end.first; iter != begin_end.second; ++iter)
285  {
286  RawTower *tmpt = towers->getTower(iter->first);
287  double e = tmpt->get_energy();
288  RawTowerGeom *tgeo =
289  towergeom->get_tower_geometry(tmpt->get_id());
290  double phi = tgeo->get_phi();
291  if (phi > M_PI) phi = phi - 2. * M_PI; // correct the cluster phi for slat geometry which is 0-2pi (L. Xue)
292  if (phi < 0.)
293  {
294  phi = phi + 2. * M_PI; // shift phi range for correct mean calculation
295  }
296  mean += e * phi;
297  }
298  mean = mean / sum;
299  if (mean > M_PI)
300  {
301  mean = mean - 2. * M_PI; // shift back
302  }
303 
304  cluster->set_phi(mean);
305 
306  return true; // mean phi was corrected
307 }
308 
310 {
312 }
313 
315 {
316  PHNodeIterator iter(topNode);
317 
318  // Grab the cEMC node
319  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
320  if (!dstNode)
321  {
322  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
323  throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
324  }
325 
326  PHNodeIterator dstiter(dstNode);
327  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
328  if (!DetNode)
329  {
330  DetNode = new PHCompositeNode(detector);
331  dstNode->addNode(DetNode);
332  }
333 
335  ClusterNodeName = "CLUSTER_" + detector;
337  DetNode->addNode(clusterNode);
338 }