ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4ZDCRawTowerBuilderByHitIndex.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4ZDCRawTowerBuilderByHitIndex.cc
2 #include "EICG4ZDCdetid.h"
3 
6 
7 #include <eiczdcbase/RawTowerZDCGeom.h> // for RawTowerGeom
9 #include <eiczdcbase/RawTowerZDCGeomContainer.h> // for RawTowerGeomContainer
10 
11 #include <g4main/PHG4Hit.h>
13 
15 #include <fun4all/SubsysReco.h> // for SubsysReco
16 
17 #include <phool/PHCompositeNode.h>
18 #include <phool/PHIODataNode.h>
19 #include <phool/PHNode.h> // for PHNode
20 #include <phool/PHNodeIterator.h>
21 #include <phool/PHObject.h> // for PHObject
22 #include <phool/getClass.h>
23 #include <phool/phool.h> // for PHWHERE
24 
25 #include <TRotation.h>
26 #include <TVector3.h>
27 
28 #include <cstdlib> // for exit
29 #include <exception> // for exception
30 #include <fstream>
31 #include <iostream>
32 #include <map>
33 #include <sstream>
34 #include <stdexcept>
35 #include <utility> // for pair, make_pair
36 #include <bitset>
37 
38 using namespace std;
39 
41  : SubsysReco(name)
42  , m_Towers(nullptr)
43  , m_Geoms(nullptr)
44  , m_Detector("NONE")
45  , m_SubDetector("NONE")
46  , m_MappingTowerFile("default.txt")
47  , m_CaloId(RawTowerZDCDefs::NONE)
48  , m_Emin(1e-9)
49  , m_TowerDepth(0)
50  , m_ThicknessAbsorber(0)
51  , m_ThicknessScintilator(0)
52 {
53 }
54 
56 {
57  PHNodeIterator iter(topNode);
58 
59  // Looking for the DST node
60  PHCompositeNode *dstNode;
61  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
62  if (!dstNode)
63  {
64  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
65  exit(1);
66  }
67 
68  try
69  {
70  CreateNodes(topNode);
71  }
72  catch (std::exception &e)
73  {
74  std::cout << e.what() << std::endl;
75  //exit(1);
76  }
77 
78  try
79  {
81  }
82  catch (std::exception &e)
83  {
84  std::cout << e.what() << std::endl;
85  //exit(1);
86  }
87 
89 }
90 
92 {
93  // get hits
94  string NodeNameHits = "G4HIT_" + m_Detector;
95  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, NodeNameHits);
96  if (!g4hit)
97  {
98  cout << "Could not locate g4 hit node " << NodeNameHits << endl;
99  exit(1);
100  }
101 
102  // loop over all hits in the event
104  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
105 
106  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
107  {
108  PHG4Hit *g4hit_i = hiter->second;
109 
110  // Don't include hits with zero energy
111  if (g4hit_i->get_edep() <= 0 && g4hit_i->get_edep() != -1) continue;
112  if (g4hit_i->get_hit_type() != m_SubDetID) continue;
113 
114  /* encode CaloTowerID from i (xid), j (yid) index of tower / hit and calorimeter ID */
115  int layerid = g4hit_i->get_layer();
116  int towerid = -1;
118  for(auto it = m_TowerIDtoLayerIDMap.begin(); it!=m_TowerIDtoLayerIDMap.end();it++)
119  if(layerid >= it->second) towerid = it->first;
120  }else {
121  for(auto it = m_TowerIDtoLayerIDMap.begin(); it!=m_TowerIDtoLayerIDMap.end(); it++){
122  if(layerid == it->second){
123  towerid = it->first;
124  break;
125  }
126  }
127  }
128 
130  g4hit_i->get_index_i(),
131  g4hit_i->get_index_j(),
132  towerid
133  );
134 
135 
136  /* add the energy to the corresponding tower */
137  RawTowerZDCv1 *tower = dynamic_cast<RawTowerZDCv1 *>(m_Towers->getTower(calotowerid));
138  if (!tower)
139  {
140  tower = new RawTowerZDCv1(calotowerid);
141  tower->set_energy(0);
142  m_Towers->AddTower(tower->get_id(), tower);
143  if (Verbosity() > 2)
144  {
145  cout<<m_SubDetector<<" : L = "<<layerid<<" , T = "<<towerid<<endl;
146 
147  std::cout << "in: " << g4hit_i->get_index_i() << "\t" << g4hit_i->get_index_j() << "\t" << towerid << std::endl;
148  std::cout << "decoded: " << tower->get_bineta() << "\t" << tower->get_binphi() << "\t" << tower->get_binl() << std::endl;
149  }
150  }
151 
152  double hit_energy = 0;
154  else hit_energy = g4hit_i->get_edep();
155 
156  tower->add_ecell(layerid, hit_energy);
157  tower->set_energy(tower->get_energy() + hit_energy);
158  tower->add_eshower(g4hit_i->get_shower_id(), hit_energy);
159 
160  }
161 
162  float towerE = 0.;
163 
164  if (Verbosity())
165  {
166  towerE = m_Towers->getTotalEdep();
167  }
168 
170  if (Verbosity())
171  {
172  std::cout << "storing towers: "<< m_Towers->size() << std::endl;
173  cout << "Energy lost by dropping towers with less than " << m_Emin
174  << " energy, lost energy: " << towerE - m_Towers->getTotalEdep() << endl;
175  m_Towers->identify();
176  std::cout<<"Total energy: "<<m_Towers->getTotalEdep()<<std::endl;
179  for (iter = begin_end.first; iter != begin_end.second; ++iter)
180  {
181  iter->second->identify();
182  }
183  }
184 
186 }
187 
189 {
191 }
192 
194 {
195  m_Detector = d;
196 
197 }
198 
200 {
201 
202  m_SubDetector = d;
203 
204  if(d == "ZDC_Crystal")
206  else if(d == "ZDC_SiPixel")
208  else if(d == "ZDC_SiPad")
210  else if(d == "ZDC_Sci")
212  else
213  {
214  std::cout<<"Invalid ZDC subdetector name in EICG4ZDCRawTowerBuilderByHitIndex"
215  <<std::endl;
216  exit(1);
217  }
218 
220 
221 }
222 
224 {
225  PHNodeIterator iter(topNode);
226  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
227  if (!runNode)
228  {
229  std::cerr << PHWHERE << "Run Node missing, doing nothing." << std::endl;
230  throw std::runtime_error("Failed to find Run node in EICG4ZDCRawTowerBuilderByHitIndex::CreateNodes");
231  }
232 
233  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
234  if (!dstNode)
235  {
236  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
237  throw std::runtime_error("Failed to find DST node in EICG4ZDCRawTowerBuilderByHitIndex::CreateNodes");
238  }
239 
240  // Create the tower geometry node on the tree
242  string NodeNameTowerGeometries = "TOWERGEOM_" + m_SubDetector;
243 
244  PHIODataNode<PHObject> *geomNode = new PHIODataNode<PHObject>(m_Geoms, NodeNameTowerGeometries, "PHObject");
245  runNode->addNode(geomNode);
246 
247  // Find detector node (or create new one if not found)
248  PHNodeIterator dstiter(dstNode);
249  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
250  "PHCompositeNode", m_SubDetector));
251  if (!DetNode)
252  {
253  DetNode = new PHCompositeNode(m_SubDetector);
254  dstNode->addNode(DetNode);
255  }
256 
257  // Create the tower nodes on the tree
259  string NodeNameTowers;
260  if (m_SimTowerNodePrefix.empty())
261  {
262  // no prefix, consistent with older convension
263  NodeNameTowers = "TOWER_" + m_SubDetector;
264  }
265  else
266  {
267  NodeNameTowers = "TOWER_" + m_SimTowerNodePrefix + "_" + m_SubDetector;
268  }
269 
270  PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_Towers, NodeNameTowers, "PHObject");
271  DetNode->addNode(towerNode);
272 
273  return;
274 }
275 
277 {
278  /* Stream to read table from file */
279  ifstream istream_mapping;
280 
281  if(m_MappingTowerFile.find(m_SubDetector)==std::string::npos){
282  cerr<<"ZDCTowerGeomManager::ReadGeometryFileFromTable - ERROR unmatched mapping file: SubDetector "<<m_SubDetector<<" : mapping file : "<<m_MappingTowerFile<<endl;
283  }
284 
285  /* Open the datafile, if it won't open return an error */
286  if (!istream_mapping.is_open())
287  {
288  istream_mapping.open(m_MappingTowerFile);
289  if (!istream_mapping)
290  {
291  cerr << "ZDCTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << m_MappingTowerFile << endl;
292  exit(1);
293  }
294  }
295 
296  string line_mapping;
297 
298  double RotAroundX=0.;
299  double RotAroundY=0.;
300  double RotAroundZ=0.;
301  double GlobalPlaceInX=0.;
302  double GlobalPlaceInY=0.;
303  double GlobalPlaceInZ=0.;
304 
305  m_TowerIDtoLayerIDMap.clear();
306 
307  while (getline(istream_mapping, line_mapping))
308  {
309  /* Skip lines starting with / including a '#' */
310  if (line_mapping.find("#") != string::npos)
311  {
312  if (Verbosity() > 0)
313  {
314  cout << "EICG4ZDCRawTowerBuilderByHitIndex: SKIPPING line in mapping file: " << line_mapping << endl;
315  }
316  continue;
317  }
318  istringstream iss(line_mapping);
319 
320  /* If line starts with keyword Tower, add to tower positions */
321  if (line_mapping.find("Tower ") != string::npos)
322  {
323  unsigned idx_T, idx_x, idx_y;
324  double pos_x, pos_y, pos_z;
325  double size_x, size_y, size_z;
326  string dummys;
327 
328  /* read string- break if error */
329  if (!(iss >> dummys >> idx_T >> idx_x >> idx_y >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z))
330  {
331  cerr << "ERROR in EICG4ZDCRawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << endl;
332  exit(1);
333  }
334 
335  /* Construct unique Tower ID */
336  unsigned int temp_id = RawTowerZDCDefs::encode_towerid_zdc(m_CaloId, idx_x, idx_y, idx_T);
337 
338  /* Create tower geometry object */
339  RawTowerZDCGeomv1 *temp_geo = new RawTowerZDCGeomv1(temp_id);
340  temp_geo->set_center_x(pos_x);
341  temp_geo->set_center_y(pos_y);
342  temp_geo->set_center_z(pos_z);
343  temp_geo->set_size_x(size_x);
344  temp_geo->set_size_y(size_y);
345  temp_geo->set_size_z(size_z);
346 
347  /* Insert this tower into position map */
348  m_Geoms->add_tower_geometry(temp_geo);
349 
350  }
351 
352  /* If line does not start with keyword Tower, read as parameter */
353  else if (line_mapping.find("iT_iL ")!=string::npos)
354  {
355  int towerid, layerid;
356  string dummys;
357  if(!(iss >> dummys >> towerid >> layerid))
358  {
359  cerr<<"ERROR in EICG4ZDCRawTowerBuilderByHitIndex: Failed to read line in mapping file "<< m_MappingTowerFile <<endl;
360  exit(1);
361  }
362  m_TowerIDtoLayerIDMap.insert(make_pair(towerid,layerid));
363  }
364  else
365  {
366  /* If this line is not a comment and not a tower, save parameter as string / value. */
367  string parname;
368  double parval;
369 
370  /* read string- break if error */
371  if (!(iss >> parname >> parval))
372  {
373  cerr << "ERROR in EICG4ZDCRawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << endl;
374  exit(1);
375  }
376 
377  m_GlobalParameterMap.insert(make_pair(parname, parval));
378 
379  /* Update member variables for global parameters based on parsed parameter file */
380  std::map<string, double>::iterator parit;
381 
382  parit = m_GlobalParameterMap.find("Gx0");
383  if (parit != m_GlobalParameterMap.end())
384  GlobalPlaceInX = parit->second;
385 
386  parit = m_GlobalParameterMap.find("Gy0");
387  if (parit != m_GlobalParameterMap.end())
388  GlobalPlaceInY = parit->second;
389 
390  parit = m_GlobalParameterMap.find("Gz0");
391  if (parit != m_GlobalParameterMap.end())
392  GlobalPlaceInZ = parit->second;
393 
394  parit = m_GlobalParameterMap.find("Grot_x");
395  if (parit != m_GlobalParameterMap.end())
396  RotAroundX = parit->second;
397 
398  parit = m_GlobalParameterMap.find("Grot_y");
399  if (parit != m_GlobalParameterMap.end())
400  RotAroundY = parit->second;
401 
402  parit = m_GlobalParameterMap.find("Grot_z");
403  if (parit != m_GlobalParameterMap.end())
404  RotAroundZ = parit->second;
405 
406  }
407  }
408 
409  /* Correct tower geometries for global calorimter translation / rotation
410  * after reading parameters from file */
412 
413  for (RawTowerZDCGeomContainer::ConstIterator it = all_towers.first;
414  it != all_towers.second; ++it)
415  {
416  double x_temp = it->second->get_center_x();
417  double y_temp = it->second->get_center_y();
418  double z_temp = it->second->get_center_z();
419 
420  TVector3 v_temp_r(x_temp, y_temp, z_temp);
421 
422  /* Rotation */
423  v_temp_r.RotateX(RotAroundX);
424  v_temp_r.RotateY(RotAroundY);
425  v_temp_r.RotateZ(RotAroundZ);
426 
427  /* Translation */
428  double x_temp_rt = v_temp_r.X() + GlobalPlaceInX;
429  double y_temp_rt = v_temp_r.Y() + GlobalPlaceInY;
430  double z_temp_rt = v_temp_r.Z() + GlobalPlaceInZ;
431 
432  /* Update tower geometry object */
433  it->second->set_center_x(x_temp_rt);
434  it->second->set_center_y(y_temp_rt);
435  it->second->set_center_z(z_temp_rt);
436 
437  if (Verbosity() > 2)
438  {
439  cout << "* Local tower x y z : " << x_temp << " " << y_temp << " " << z_temp << endl;
440  cout << "* Globl tower x y z : " << x_temp_rt << " " << y_temp_rt << " " << z_temp_rt << endl;
441  }
442  }
443 
444  return true;
445 }