Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawClusterBuilderGraph.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawClusterBuilderGraph.cc
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/PHCompositeNode.h>
19 #include <phool/PHIODataNode.h>
20 #include <phool/PHNode.h>
21 #include <phool/PHNodeIterator.h>
22 #include <phool/PHObject.h>
23 #include <phool/getClass.h>
24 #include <phool/phool.h>
25 
26 #include <algorithm> // for sort
27 #include <cassert>
28 #include <cmath>
29 #include <exception>
30 #include <iostream>
31 #include <map>
32 #include <stdexcept>
33 #include <utility>
34 #include <vector>
35 
36 // this is just a helper class which enables us to handle rollovers
37 // when checking for adjacent towers, it requires one bit of
38 // information (the total number of phibins) which
39 // is not in the tower class
40 class twrs
41 {
42  public:
43  explicit twrs(RawTower * /*rt*/);
44  virtual ~twrs() = default;
45  bool is_adjacent(const twrs & /*tower*/);
46  void set_id(const int i)
47  {
48  id = i;
49  }
50  int get_id() const
51  {
52  return id;
53  }
54  void set_maxphibin(const int i)
55  {
56  maxphibin = i;
57  }
58  int get_maxphibin() const
59  {
60  return maxphibin;
61  }
62  int get_bineta() const
63  {
64  return bineta;
65  }
66  int get_binphi() const
67  {
68  return binphi;
69  }
70 
71  private:
72  int bineta;
73  int binphi;
74  int maxphibin;
76 };
77 
79  : maxphibin(-10)
80  , id(-1)
81 {
82  bineta = rt->get_bineta();
83  binphi = rt->get_binphi();
84 }
85 
87 {
88  if (bineta - 1 <= tower.get_bineta() && tower.get_bineta() <= bineta + 1)
89  {
90  if (binphi - 1 <= tower.get_binphi() && tower.get_binphi() <= binphi + 1)
91  {
92  return true;
93  }
94  // cluster through the phi-wraparound
95  else if (((tower.get_binphi() == maxphibin - 1) && (binphi == 0)) ||
96  ((tower.get_binphi() == 0) && (binphi == maxphibin - 1)))
97  {
98  return true;
99  }
100  }
101 
102  return false;
103 }
104 
105 bool operator<(const twrs &a, const twrs &b)
106 {
107  if (a.get_bineta() != b.get_bineta())
108  {
109  return a.get_bineta() < b.get_bineta();
110  }
111  return a.get_binphi() < b.get_binphi();
112 }
113 
115  : SubsysReco(name)
116  , _clusters(nullptr)
117  , _min_tower_e(0.0)
118  , chkenergyconservation(0)
119  , detector("NONE")
120 {
121 }
122 
124 {
125  try
126  {
127  CreateNodes(topNode);
128  }
129  catch (std::exception &e)
130  {
131  std::cout << PHWHERE << ": " << e.what() << std::endl;
132  throw;
133  }
134 
136 }
137 
139 {
140  std::string towernodename = "TOWER_CALIB_" + detector;
141  // Grab the towers
142  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, towernodename);
143  if (!towers)
144  {
145  std::cout << PHWHERE << ": Could not find node " << towernodename << std::endl;
147  }
148  std::string towergeomnodename = "TOWERGEOM_" + detector;
149  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
150  if (!towergeom)
151  {
152  std::cout << PHWHERE << ": Could not find node " << towergeomnodename << std::endl;
154  }
155  // make the list of towers above threshold
156  std::vector<twrs> towerVector;
157  RawTowerContainer::ConstRange begin_end = towers->getTowers();
158  RawTowerContainer::ConstIterator itr = begin_end.first;
159  for (; itr != begin_end.second; ++itr)
160  {
161  RawTower *tower = itr->second;
162  RawTowerDefs::keytype towerid = itr->first;
163  if (tower->get_energy() > _min_tower_e)
164  {
165  twrs twr(tower);
166  twr.set_maxphibin(towergeom->get_phibins());
167  twr.set_id(towerid);
168  towerVector.push_back(twr);
169  }
170  }
171 
172  // cluster the towers
173  std::multimap<int, twrs> clusteredTowers;
174  PHMakeGroups(towerVector, clusteredTowers);
175 
176  RawCluster *cluster = nullptr;
177  int last_id = -1;
178  std::multimap<int, twrs>::iterator ctitr = clusteredTowers.begin();
179  std::multimap<int, twrs>::iterator lastct = clusteredTowers.end();
180  for (; ctitr != lastct; ++ctitr)
181  {
182  int clusterid = ctitr->first;
183 
184  if (last_id != clusterid)
185  {
186  // new cluster
187  cluster = new RawClusterv1();
188  _clusters->AddCluster(cluster);
189 
190  last_id = clusterid;
191  }
192  assert(cluster);
193 
194  const twrs &tmptower = ctitr->second;
195  RawTower *rawtower = towers->getTower(tmptower.get_id());
196 
197  const double e = rawtower->get_energy();
198  cluster->addTower(rawtower->get_id(), e);
199  }
200 
201  for (const auto &cluster_pair : _clusters->getClustersMap())
202  {
203  RawClusterDefs::keytype clusterid = cluster_pair.first;
204  RawCluster *clusterA = cluster_pair.second;
205 
206  assert(clusterA);
207  assert(clusterA->get_id() == clusterid);
208 
209  double sum_x(0);
210  double sum_y(0);
211  double sum_z(0);
212  double sum_e(0);
213 
214  for (const auto tower_pair : clusterA->get_towermap())
215  {
216  const RawTower *rawtower = towers->getTower(tower_pair.first);
217  const RawTowerGeom *rawtowergeom = towergeom->get_tower_geometry(tower_pair.first);
218 
219  assert(rawtower);
220  assert(rawtowergeom);
221  const double e = rawtower->get_energy();
222 
223  sum_e += e;
224 
225  if (e > 0)
226  {
227  sum_x += e * rawtowergeom->get_center_x();
228  sum_y += e * rawtowergeom->get_center_y();
229  sum_z += e * rawtowergeom->get_center_z();
230  }
231  } // for (const auto tower_pair : clusterA->get_towermap())
232 
233  clusterA->set_energy(sum_e);
234 
235  if (sum_e > 0)
236  {
237  sum_x /= sum_e;
238  sum_y /= sum_e;
239  sum_z /= sum_e;
240 
241  clusterA->set_r(sqrt(sum_y * sum_y + sum_x * sum_x));
242  clusterA->set_phi(atan2(sum_y, sum_x));
243  clusterA->set_z(sum_z);
244  }
245 
246  if (Verbosity() > 1)
247  {
248  std::cout << "RawClusterBuilderGraph constucted ";
249  clusterA->identify();
250  }
251  } // for (const auto & cluster_pair : _clusters->getClustersMap())
252 
254  {
255  double ecluster = _clusters->getTotalEdep();
256  double etower = towers->getTotalEdep();
257  if (ecluster > 0)
258  {
259  if (fabs(etower - ecluster) / ecluster > 1e-9)
260  {
261  std::cout << "energy conservation violation: ETower: " << etower
262  << " ECluster: " << ecluster
263  << " diff: " << etower - ecluster << std::endl;
264  }
265  }
266  else
267  {
268  if (etower != 0)
269  {
270  std::cout << "energy conservation violation: ETower: " << etower
271  << " ECluster: " << ecluster << std::endl;
272  }
273  }
274  }
276 }
277 
279 {
281 }
282 
284 {
285  PHNodeIterator iter(topNode);
286 
287  // Grab the cEMC node
288  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
289  if (!dstNode)
290  {
291  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
292  throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
293  }
294 
295  PHNodeIterator dstiter(dstNode);
296  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
297  if (!DetNode)
298  {
299  DetNode = new PHCompositeNode(detector);
300  dstNode->addNode(DetNode);
301  }
302 
304  ClusterNodeName = "CLUSTER_" + detector;
306  DetNode->addNode(clusterNode);
307 }