ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerBuilderByHitIndexBECAL.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerBuilderByHitIndexBECAL.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/RawTowerGeomContainer_Cylinderv1.h>
11 #include <calobase/RawTowerGeomv4.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_Emin(1e-6)
49 {
50 }
51 
53 {
54  PHNodeIterator iter(topNode);
55 
56  // Looking for the DST node
57  PHCompositeNode *dstNode;
58  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
59  if (!dstNode)
60  {
61  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
62  exit(1);
63  }
64 
65  try
66  {
67  CreateNodes(topNode);
68  }
69  catch (std::exception &e)
70  {
71  std::cout << e.what() << std::endl;
72  //exit(1);
73  }
74 
75  try
76  {
78  }
79  catch (std::exception &e)
80  {
81  std::cout << e.what() << std::endl;
82  //exit(1);
83  }
84 
86 }
87 
89 {
90  // get hits
91  string NodeNameHits = "G4HIT_" + m_Detector;
92  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, NodeNameHits);
93  if (!g4hit)
94  {
95  cout << "Could not locate g4 hit node " << NodeNameHits << endl;
96  exit(1);
97  }
98 
99  // loop over all hits in the event
101  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
102 
103  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
104  {
105  PHG4Hit *g4hit_i = hiter->second;
106 
107  // Don't include hits with zero energy
108  if (g4hit_i->get_edep() <= 0 && g4hit_i->get_edep() != -1) continue;
109 
110  /* encode CaloTowerID from j, k indexBECAL of tower / hit and calorimeter ID */
112  g4hit_i->get_index_j(),
113  g4hit_i->get_index_k());
114 
115  /* add the energy to the corresponding tower */
116  RawTowerv1 *tower = dynamic_cast<RawTowerv1 *>(m_Towers->getTower(calotowerid));
117  if (!tower)
118  {
119  tower = new RawTowerv1(calotowerid);
120  tower->set_energy(0);
121  m_Towers->AddTower(tower->get_id(), tower);
122  }
123 
124  tower->add_ecell((g4hit_i->get_index_j() << 16) + g4hit_i->get_index_k(), g4hit_i->get_light_yield());
125  tower->set_energy(tower->get_energy() + g4hit_i->get_light_yield());
126  tower->add_eshower(g4hit_i->get_shower_id(), g4hit_i->get_edep());
127  }
128 
129  float towerE = 0.;
130 
131  if (Verbosity())
132  {
133  towerE = m_Towers->getTotalEdep();
134  std::cout << "towers before compression: " << m_Towers->size() << "\t" << m_Detector << std::endl;
135  }
137  if (Verbosity())
138  {
139  std::cout << "storing towers: " << m_Towers->size() << std::endl;
140  cout << "Energy lost by dropping towers with less than " << m_Emin
141  << " energy, lost energy: " << towerE - m_Towers->getTotalEdep() << endl;
142  m_Towers->identify();
145  for (iter = begin_end.first; iter != begin_end.second; ++iter)
146  {
147  iter->second->identify();
148  }
149  }
150 
152 }
153 
155 {
157 }
158 
160 {
161  m_Detector = d;
163 }
164 
166 {
167  PHNodeIterator iter(topNode);
168  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
169  if (!runNode)
170  {
171  std::cerr << PHWHERE << "Run Node missing, doing nothing." << std::endl;
172  throw std::runtime_error("Failed to find Run node in RawTowerBuilderByHitIndexBECAL::CreateNodes");
173  }
174 
175  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
176  if (!dstNode)
177  {
178  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
179  throw std::runtime_error("Failed to find DST node in RawTowerBuilderByHitIndexBECAL::CreateNodes");
180  }
181 
182  // Create the tower geometry node on the tree
184 
185  string NodeNameTowerGeometries = "TOWERGEOM_" + m_Detector;
186 
187  PHIODataNode<PHObject> *geomNode = new PHIODataNode<PHObject>(m_Geoms, NodeNameTowerGeometries, "PHObject");
188  runNode->addNode(geomNode);
189 
190  // Find detector node (or create new one if not found)
191  PHNodeIterator dstiter(dstNode);
192  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
193  "PHCompositeNode", m_Detector));
194  if (!DetNode)
195  {
196  DetNode = new PHCompositeNode(m_Detector);
197  dstNode->addNode(DetNode);
198  }
199 
200  // Create the tower nodes on the tree
202  string NodeNameTowers;
203  if (m_SimTowerNodePrefix.empty())
204  {
205  // no prefix, consistent with older convension
206  NodeNameTowers = "TOWER_" + m_Detector;
207  }
208  else
209  {
210  NodeNameTowers = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
211  }
212 
213  PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_Towers, NodeNameTowers, "PHObject");
214  DetNode->addNode(towerNode);
215 
216  return;
217 }
218 
220 {
221  /* Stream to read table from file */
222  ifstream istream_mapping;
223 
224  /* Open the datafile, if it won't open return an error */
225  if (!istream_mapping.is_open())
226  {
227  istream_mapping.open(m_MappingTowerFile);
228  if (!istream_mapping)
229  {
230  cerr << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << m_MappingTowerFile << endl;
231  exit(1);
232  }
233  }
234 
235  string line_mapping;
236 
237  while (getline(istream_mapping, line_mapping))
238  {
239  /* Skip lines starting with / including a '#' */
240  if (line_mapping.find("#") != string::npos)
241  {
242  if (Verbosity() > 0)
243  {
244  cout << "RawTowerBuilderByHitIndexBECAL: SKIPPING line in mapping file: " << line_mapping << endl;
245  }
246  continue;
247  }
248 
249  std::istringstream iss(line_mapping);
250 
251  if (line_mapping.find("BECALtower ") != string::npos)
252  {
253  unsigned idphi_j, ideta_k;
254  double cx, cy, cz;
255  double rot_z, rot_y, rot_x;
256  double size_z, size_y1, size_x1, size_y2, size_x2, p_Theta;
257  std::string dummys;
258 
259  if (!(iss >> dummys >> ideta_k >> idphi_j >> size_x1 >> size_x2 >> size_y1 >> size_y2 >> size_z >> p_Theta >> cx >> cy >> cz >> rot_x >> rot_y >> rot_z))
260  {
261  std::cout << "ERROR in PHG4ForwardHcalDetector: Failed to read line in mapping file " << m_MappingTowerFile << std::endl;
262  exit(1);
263  }
264 
265  /* Construct unique Tower ID */
266  unsigned int temp_id = RawTowerDefs::encode_towerid(m_CaloId, ideta_k, idphi_j);
267 
268  /* Create tower geometry object */
269  RawTowerGeom *temp_geo = new RawTowerGeomv4(temp_id);
270  temp_geo->set_center_x(cx);
271  temp_geo->set_center_y(cy);
272  temp_geo->set_center_z(cz);
273  temp_geo->set_roty(rot_y);
274  temp_geo->set_rotz(rot_z);
275 
276  m_Geoms->add_tower_geometry(temp_geo);
277  }
278  else
279  {
280  string parname;
281  double parval;
282  if (!(iss >> parname >> parval))
283  {
284  cout << "ERROR in PHG4BarrelCalorimeterDetector: Failed to read line in mapping file " << endl;
285  }
286 
287  m_GlobalParameterMap.insert(make_pair(parname, parval));
288 
289  std::map<string, double>::iterator parit;
290 
291  parit = m_GlobalParameterMap.find("thickness_wall");
292  if (parit != m_GlobalParameterMap.end())
293  {
294  thickness_wall = parit->second; // in cm
295  }
296 
297  parit = m_GlobalParameterMap.find("radius");
298  if (parit != m_GlobalParameterMap.end())
299  {
300  radius = parit->second; // in cm
301  }
302 
303  parit = m_GlobalParameterMap.find("tower_length");
304  if (parit != m_GlobalParameterMap.end())
305  {
306  tower_length = parit->second; // in cm
307  }
308  }
309  }
310 
312  //cout << PHWHERE << "BECAL radius set to " << m_Geoms->get_radius() << " cm" << endl;
313 
315 
316  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
317  it != all_towers.second; ++it)
318  {
319  double x_temp = it->second->get_center_x();
320  double y_temp = it->second->get_center_y();
321  double z_temp = it->second->get_center_z();
322  double roty = it->second->get_roty();
323  double rotz = it->second->get_rotz();
324 
325  double x_tempfinal = x_temp + thickness_wall / 2 * abs(cos(roty - M_PI_2)) * cos(rotz);
326  double y_tempfinal = y_temp + thickness_wall / 2 * abs(cos(roty - M_PI_2)) * sin(rotz);
327  double z_tempfinal = z_temp + thickness_wall * sin(roty - M_PI_2);
328 
329  it->second->set_center_x(x_tempfinal);
330  it->second->set_center_y(y_tempfinal);
331  it->second->set_center_z(z_tempfinal);
332 
333  if (Verbosity() > 2)
334  {
335  cout << "*** Tower x y z : " << x_temp << " " << y_temp << " " << z_temp << endl;
336  }
337  }
338  if (Verbosity())
339  {
340  cout << "size tower geom container:" << m_Geoms->size() << "\t" << m_Detector << endl;
341  }
342 
343  return true;
344 }