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