ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerBuilderByHitIndex.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerBuilderByHitIndex.cc
2 
3 #include <calobase/RawTowerContainer.h>
4 #include <calobase/RawTowerv1.h>
5 
6 #include <calobase/RawTower.h> // for RawTower
7 #include <calobase/RawTowerDefs.h> // for convert_name_to_caloid
8 #include <calobase/RawTowerGeom.h> // for RawTowerGeom
9 #include <calobase/RawTowerGeomv3.h>
10 #include <calobase/RawTowerGeomContainer.h> // for RawTowerGeomContainer
11 #include <calobase/RawTowerGeomContainerv1.h>
12 
13 #include <g4main/PHG4Hit.h>
15 
17 #include <fun4all/SubsysReco.h> // for SubsysReco
18 
19 #include <phool/PHCompositeNode.h>
20 #include <phool/PHIODataNode.h>
21 #include <phool/PHNode.h> // for PHNode
22 #include <phool/PHNodeIterator.h>
23 #include <phool/PHObject.h> // for PHObject
24 #include <phool/getClass.h>
25 #include <phool/phool.h> // for PHWHERE
26 
27 #include <TRotation.h>
28 #include <TVector3.h>
29 
30 #include <cstdlib> // for exit
31 #include <exception> // for exception
32 #include <fstream>
33 #include <iostream>
34 #include <map>
35 #include <sstream>
36 #include <stdexcept>
37 #include <utility> // for pair, make_pair
38 
39 using namespace std;
40 
42  : SubsysReco(name)
43  , m_Towers(nullptr)
44  , m_Geoms(nullptr)
45  , m_Detector("NONE")
46  , m_MappingTowerFile("default.txt")
47  , m_CaloId(RawTowerDefs::NONE)
48  , m_GlobalPlaceInX(0)
49  , m_GlobalPlaceInY(0)
50  , m_GlobalPlaceInZ(0)
51  , m_RotInX(0)
52  , m_RotInY(0)
53  , m_RotInZ(0)
54  , m_Emin(1e-6)
55 {
56 }
57 
59 {
60  PHNodeIterator iter(topNode);
61 
62  // Looking for the DST node
63  PHCompositeNode *dstNode;
64  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
65  if (!dstNode)
66  {
67  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
68  exit(1);
69  }
70 
71  try
72  {
73  CreateNodes(topNode);
74  }
75  catch (std::exception &e)
76  {
77  std::cout << e.what() << std::endl;
78  //exit(1);
79  }
80 
81  try
82  {
84  }
85  catch (std::exception &e)
86  {
87  std::cout << e.what() << std::endl;
88  //exit(1);
89  }
90 
92 }
93 
95 {
96  // get hits
97  string NodeNameHits = "G4HIT_" + m_Detector;
98  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, NodeNameHits);
99  if (!g4hit)
100  {
101  cout << "Could not locate g4 hit node " << NodeNameHits << endl;
102  exit(1);
103  }
104 
105  // loop over all hits in the event
107  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
108 
109  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
110  {
111  PHG4Hit *g4hit_i = hiter->second;
112 
113  // Don't include hits with zero energy
114  if (g4hit_i->get_edep() <= 0 && g4hit_i->get_edep() != -1) continue;
115 
116  /* encode CaloTowerID from j, k index of tower / hit and calorimeter ID */
118  g4hit_i->get_index_j(),
119  g4hit_i->get_index_k());
120 
121  /* add the energy to the corresponding tower */
122  RawTowerv1 *tower = dynamic_cast<RawTowerv1 *>(m_Towers->getTower(calotowerid));
123  if (!tower)
124  {
125  tower = new RawTowerv1(calotowerid);
126  tower->set_energy(0);
127  m_Towers->AddTower(tower->get_id(), tower);
128  }
129 
130  tower->add_ecell((g4hit_i->get_index_j() << 16) + g4hit_i->get_index_k(), g4hit_i->get_light_yield());
131  tower->set_energy(tower->get_energy() + g4hit_i->get_light_yield());
132  tower->add_eshower(g4hit_i->get_shower_id(), g4hit_i->get_edep());
133  }
134 
135  float towerE = 0.;
136 
137  if (Verbosity())
138  {
139  towerE = m_Towers->getTotalEdep();
140  std::cout << "towers before compression: "<< m_Towers->size() << "\t" << m_Detector << std::endl;
141  }
143  if (Verbosity())
144  {
145  std::cout << "storing towers: "<< m_Towers->size() << std::endl;
146  cout << "Energy lost by dropping towers with less than " << m_Emin
147  << " energy, lost energy: " << towerE - m_Towers->getTotalEdep() << endl;
148  m_Towers->identify();
151  for (iter = begin_end.first; iter != begin_end.second; ++iter)
152  {
153  iter->second->identify();
154  }
155  }
156 
158 }
159 
161 {
163 }
164 
165 void RawTowerBuilderByHitIndex::Detector(const std::string &d)
166 {
167  m_Detector = d;
169 }
170 
172 {
173  PHNodeIterator iter(topNode);
174  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
175  if (!runNode)
176  {
177  std::cerr << PHWHERE << "Run Node missing, doing nothing." << std::endl;
178  throw std::runtime_error("Failed to find Run node in RawTowerBuilderByHitIndex::CreateNodes");
179  }
180 
181  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
182  if (!dstNode)
183  {
184  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
185  throw std::runtime_error("Failed to find DST node in RawTowerBuilderByHitIndex::CreateNodes");
186  }
187 
188  // Create the tower geometry node on the tree
190  string NodeNameTowerGeometries = "TOWERGEOM_" + m_Detector;
191 
192  PHIODataNode<PHObject> *geomNode = new PHIODataNode<PHObject>(m_Geoms, NodeNameTowerGeometries, "PHObject");
193  runNode->addNode(geomNode);
194 
195  // Find detector node (or create new one if not found)
196  PHNodeIterator dstiter(dstNode);
197  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
198  "PHCompositeNode", m_Detector));
199  if (!DetNode)
200  {
201  DetNode = new PHCompositeNode(m_Detector);
202  dstNode->addNode(DetNode);
203  }
204 
205  // Create the tower nodes on the tree
207  string NodeNameTowers;
208  if (m_SimTowerNodePrefix.empty())
209  {
210  // no prefix, consistent with older convension
211  NodeNameTowers = "TOWER_" + m_Detector;
212  }
213  else
214  {
215  NodeNameTowers = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
216  }
217 
218  PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_Towers, NodeNameTowers, "PHObject");
219  DetNode->addNode(towerNode);
220 
221  return;
222 }
223 
225 {
226  /* Stream to read table from file */
227  ifstream istream_mapping;
228 
229  /* Open the datafile, if it won't open return an error */
230  if (!istream_mapping.is_open())
231  {
232  istream_mapping.open(m_MappingTowerFile);
233  if (!istream_mapping)
234  {
235  cerr << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << m_MappingTowerFile << endl;
236  exit(1);
237  }
238  }
239 
240  string line_mapping;
241 
242  while (getline(istream_mapping, line_mapping))
243  {
244  /* Skip lines starting with / including a '#' */
245  if (line_mapping.find("#") != string::npos)
246  {
247  if (Verbosity() > 0)
248  {
249  cout << "RawTowerBuilderByHitIndex: SKIPPING line in mapping file: " << line_mapping << endl;
250  }
251  continue;
252  }
253 
254  istringstream iss(line_mapping);
255 
256  /* If line starts with keyword Tower, add to tower positions */
257  if (line_mapping.find("Tower ") != string::npos)
258  {
259  unsigned idx_j, idx_k, idx_l;
260  double pos_x, pos_y, pos_z;
261  double size_x, size_y, size_z;
262  double rot_x, rot_y, rot_z;
263  double type;
264  string dummys;
265 
266  /* read string- break if error */
267  if (!(iss >> dummys >> type >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z >> rot_x >> rot_y >> rot_z))
268  {
269  cerr << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << endl;
270  exit(1);
271  }
272 
273  /* Construct unique Tower ID */
274  unsigned int temp_id = RawTowerDefs::encode_towerid(m_CaloId, idx_j, idx_k);
275 
276  /* Create tower geometry object */
277  RawTowerGeom *temp_geo = new RawTowerGeomv3(temp_id);
278  temp_geo->set_center_x(pos_x);
279  temp_geo->set_center_y(pos_y);
280  temp_geo->set_center_z(pos_z);
281  temp_geo->set_size_x(size_x);
282  temp_geo->set_size_y(size_y);
283  temp_geo->set_size_z(size_z);
284  temp_geo->set_tower_type((int) type);
285 
286  /* Insert this tower into position map */
287  m_Geoms->add_tower_geometry(temp_geo);
288  }
289  /* If line does not start with keyword Tower, read as parameter */
290  else
291  {
292  /* If this line is not a comment and not a tower, save parameter as string / value. */
293  string parname;
294  double parval;
295 
296  /* read string- break if error */
297  if (!(iss >> parname >> parval))
298  {
299  cerr << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << endl;
300  exit(1);
301  }
302 
303  m_GlobalParameterMap.insert(make_pair(parname, parval));
304  }
305  }
306 
307  /* Update member variables for global parameters based on parsed parameter file */
308  std::map<string, double>::iterator parit;
309 
310  parit = m_GlobalParameterMap.find("Gx0");
311  if (parit != m_GlobalParameterMap.end())
312  m_GlobalPlaceInX = parit->second;
313 
314  parit = m_GlobalParameterMap.find("Gy0");
315  if (parit != m_GlobalParameterMap.end())
316  m_GlobalPlaceInY = parit->second;
317 
318  parit = m_GlobalParameterMap.find("Gz0");
319  if (parit != m_GlobalParameterMap.end())
320  m_GlobalPlaceInZ = parit->second;
321 
322  parit = m_GlobalParameterMap.find("Grot_x");
323  if (parit != m_GlobalParameterMap.end())
324  m_RotInX = parit->second;
325 
326  parit = m_GlobalParameterMap.find("Grot_y");
327  if (parit != m_GlobalParameterMap.end())
328  m_RotInY = parit->second;
329 
330  parit = m_GlobalParameterMap.find("Grot_z");
331  if (parit != m_GlobalParameterMap.end())
332  m_RotInZ = parit->second;
333 
334  /* Correct tower geometries for global calorimter translation / rotation
335  * after reading parameters from file */
337 
338  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
339  it != all_towers.second; ++it)
340  {
341  double x_temp = it->second->get_center_x();
342  double y_temp = it->second->get_center_y();
343  double z_temp = it->second->get_center_z();
344 
345  /* Rotation */
346  TRotation rot;
347  rot.RotateX(m_RotInX);
348  rot.RotateY(m_RotInY);
349  rot.RotateZ(m_RotInZ);
350 
351  TVector3 v_temp_r(x_temp, y_temp, z_temp);
352  v_temp_r.Transform(rot);
353 
354  /* Translation */
355  double x_temp_rt = v_temp_r.X() + m_GlobalPlaceInX;
356  double y_temp_rt = v_temp_r.Y() + m_GlobalPlaceInY;
357  double z_temp_rt = v_temp_r.Z() + m_GlobalPlaceInZ;
358 
359  /* Update tower geometry object */
360  it->second->set_center_x(x_temp_rt);
361  it->second->set_center_y(y_temp_rt);
362  it->second->set_center_z(z_temp_rt);
363 
364  if (Verbosity() > 2)
365  {
366  cout << "* Local tower x y z : " << x_temp << " " << y_temp << " " << z_temp << endl;
367  cout << "* Globl tower x y z : " << x_temp_rt << " " << y_temp_rt << " " << z_temp_rt << endl;
368  }
369  }
370 
371  if (Verbosity())
372  {
373  cout << "size tower geom container:" << m_Geoms->size() << "\t" << m_Detector << endl;
374  }
375  return true;
376 }