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