ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawClusterBuilderHelper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawClusterBuilderHelper.cc
2 
3 #include <calobase/RawCluster.h>
4 #include <calobase/RawClusterContainer.h>
5 #include <calobase/RawClusterDefs.h>
6 #include <calobase/RawClusterv1.h>
7 #include <calobase/RawTower.h>
8 #include <calobase/RawTowerContainer.h>
9 #include <calobase/RawTowerDefs.h>
10 #include <calobase/RawTowerGeom.h>
11 #include <calobase/RawTowerGeomContainer.h>
12 
14 #include <fun4all/SubsysReco.h>
15 
16 #include <phool/PHCompositeNode.h>
17 #include <phool/PHIODataNode.h>
18 #include <phool/PHNode.h>
19 #include <phool/PHNodeIterator.h>
20 #include <phool/PHObject.h>
21 #include <phool/getClass.h>
22 #include <phool/phool.h>
23 
24 #include <algorithm>
25 #include <cassert>
26 #include <cmath>
27 #include <exception>
28 #include <iostream>
29 #include <map>
30 #include <stdexcept>
31 #include <string>
32 #include <utility>
33 #include <vector>
34 
36  : SubsysReco(name)
37  , _seed_e(0.5)
38  , _agg_e(0.1)
39  , detector("NONE")
40  , _clusters(nullptr)
41 {
42 }
43 
45 {
46  try
47  {
48  CreateNodes(topNode);
49  }
50  catch (std::exception &e)
51  {
52  std::cout << PHWHERE << ": " << e.what() << std::endl;
53  throw;
54  }
55 
57 }
59 {
60  std::string towernodename = "TOWER_CALIB_" + detector;
61  // Grab the towers
62  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, towernodename);
63  if (!towers)
64  {
65  std::cout << PHWHERE << ": Could not find node " << towernodename << std::endl;
67  }
68  std::string towergeomnodename = "TOWERGEOM_" + detector;
69  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
70  if (!towergeom)
71  {
72  std::cout << PHWHERE << ": Could not find node " << towergeomnodename << std::endl;
74  }
75 
76  // make the list of towers above threshold
77  std::vector<towersStrct> input_towers;
78  // towers in the current cluster
79  int towers_added = 0;
80  RawTowerContainer::ConstRange begin_end = towers->getTowers();
81  for (RawTowerContainer::ConstIterator itr = begin_end.first; itr != begin_end.second; ++itr)
82  {
83  RawTower *tower = itr->second;
84  RawTowerDefs::keytype towerid = itr->first;
85  if (tower->get_energy() > _agg_e)
86  {
87  towersStrct tempTower;
88  tempTower.tower_E = tower->get_energy();
89  tempTower.tower_iEta = tower->get_bineta();
90  tempTower.tower_iPhi = tower->get_binphi();
91  tempTower.tower_iL = 0;
92  if (towerid == RawTowerDefs::LFHCAL)
93  {
94  tempTower.tower_iL = tower->get_binl();
95  }
96  tempTower.tower_trueID = towerid; // currently unsigned -> signed, will this matter?
97  tempTower.twr = itr->second;
98  input_towers.push_back(tempTower);
99  towers_added++;
100  }
101  }
102 
103  cluster(input_towers, towers->getCalorimeterID());
104 
105  // Sum x, y, z, e
106  // from https://github.com/ECCE-EIC/coresoftware/blob/ae0526adf82f49cb8906d447411b90287de6a56e/offline/packages/CaloReco/RawClusterBuilderGraph.cc#L202
107  for (const auto &cluster_pair : _clusters->getClustersMap())
108  {
109  RawClusterDefs::keytype clusterid = cluster_pair.first;
110  RawCluster *cluster = cluster_pair.second;
111 
112  assert(cluster);
113  assert(cluster->get_id() == clusterid);
114 
115  double sum_x(0);
116  double sum_y(0);
117  double sum_z(0);
118  double sum_e(0);
119 
120  for (const auto tower_pair : cluster->get_towermap())
121  {
122  const RawTower *rawtower = towers->getTower(tower_pair.first);
123  const RawTowerGeom *rawtowergeom = towergeom->get_tower_geometry(tower_pair.first);
124 
125  assert(rawtower);
126  assert(rawtowergeom);
127  const double e = rawtower->get_energy();
128 
129  sum_e += e;
130 
131  if (e > 0)
132  {
133  sum_x += e * rawtowergeom->get_center_x();
134  sum_y += e * rawtowergeom->get_center_y();
135  sum_z += e * rawtowergeom->get_center_z();
136  }
137  } // for (const auto tower_pair : cluster->get_towermap())
138 
139  cluster->set_energy(sum_e);
140 
141  if (sum_e > 0)
142  {
143  sum_x /= sum_e;
144  sum_y /= sum_e;
145  sum_z /= sum_e;
146 
147  cluster->set_r(sqrt(sum_y * sum_y + sum_x * sum_x));
148  cluster->set_phi(atan2(sum_y, sum_x));
149 
150  cluster->set_z(sum_z);
151  }
152 
153  if (Verbosity() > 1)
154  {
155  std::cout << "RawClusterBuilderGraph constucted ";
156  cluster->identify();
157  } // for (const auto & cluster_pair : _clusters->getClustersMap())
158 
159  // The output of the cluster will be in the RawClusterContainer class
160  }
161  if (Verbosity() > 1)
162  {
163  std::cout << "Found " << _clusters->getClustersMap().size() << " clusters in " << Name() << std::endl;
164  for (const auto &cluster_pair : _clusters->getClustersMap())
165  {
166  std::cout << "\n\tnTowers: " << cluster_pair.second->getNTowers() << std::endl;
167  std::cout << "\tE: " << cluster_pair.second->get_energy() << "\tPhi: " << cluster_pair.second->get_phi() << std::endl;
168  std::cout << "\tX: " << cluster_pair.second->get_x();
169  std::cout << "\tY: " << cluster_pair.second->get_y();
170  std::cout << "\tZ: " << cluster_pair.second->get_z() << std::endl;
171  }
172  }
173 
175 }
176 
178 {
180 }
181 
183 {
184  switch (caloID)
185  {
187  return true;
188  case RawTowerDefs::FHCAL:
189  return true;
190  case RawTowerDefs::FEMC:
191  return true;
192  case RawTowerDefs::EHCAL:
193  return true;
194  case RawTowerDefs::EEMC:
195  return true;
197  return false;
199  return false;
200  case RawTowerDefs::CEMC:
201  return false;
203  return true;
205  return true;
207  return true;
208  case RawTowerDefs::BECAL:
209  return false;
210  default:
211  std::cout << "IsForwardCalorimeter: caloID " << caloID << " not defined, returning false" << std::endl;
212  return false;
213  }
214  return false;
215 }
216 
218 {
219  switch (caloID)
220  {
221  case RawTowerDefs::CEMC:
222  return 100;
224  return 64;
226  return 64;
228  return -1;
229  case RawTowerDefs::BECAL:
230  return 128;
231  default:
232  return 0;
233  }
234 }
235 
237 {
238  PHNodeIterator iter(topNode);
239 
240  // Grab the cEMC node
241  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
242  if (!dstNode)
243  {
244  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
245  throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
246  }
247 
248  PHNodeIterator dstiter(dstNode);
249  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
250  if (!DetNode)
251  {
252  DetNode = new PHCompositeNode(detector);
253  dstNode->addNode(DetNode);
254  }
255 
257  ClusterNodeName = "CLUSTER_" + detector;
259  DetNode->addNode(clusterNode);
260 }