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