ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BwdRawTowerBuilderByHitIndex.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BwdRawTowerBuilderByHitIndex.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 <pdbcalbase/PdbParameterMap.h>
26 
27 #include <TRotation.h>
28 #include <TVector3.h>
29 #include <Geant4/G4SystemOfUnits.hh>
30 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
31 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
32 #include <Geant4/G4Types.hh> // for G4double, G4int
33 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
34 #include <TSystem.h>
35 #include <Geant4/G4UnionSolid.hh>
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 
46 using namespace std;
47 
49  : SubsysReco(name)
50  , m_Towers(nullptr)
51  , m_Geoms(nullptr)
52  , m_Detector("NONE")
53  , m_MappingTowerFile("default.txt")
54  , m_CaloId(RawTowerDefs::NONE)
55  , m_GlobalPlaceInX(0)
56  , m_GlobalPlaceInY(0)
57  , m_GlobalPlaceInZ(0)
58  , m_RotInX(0)
59  , m_RotInY(0)
60  , m_RotInZ(0)
61  , m_Emin(1e-6)
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  if (g4hit_i->get_index_j() < 0 || g4hit_i->get_index_k() <0) continue;
123 
124  /* encode CaloTowerID from j, k index of tower / hit and calorimeter ID */
125 //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;
127  g4hit_i->get_index_j(),
128  g4hit_i->get_index_k());
129  /* add the energy to the corresponding tower */
130  RawTowerv1 *tower = dynamic_cast<RawTowerv1 *>(m_Towers->getTower(calotowerid));
131  if (!tower)
132  {
133  tower = new RawTowerv1(calotowerid);
134  tower->set_energy(0);
135  m_Towers->AddTower(tower->get_id(), tower);
136  }
137 
138  tower->add_ecell((g4hit_i->get_index_j() << 16) + g4hit_i->get_index_k(), g4hit_i->get_light_yield());
139  tower->set_energy(tower->get_energy() + g4hit_i->get_light_yield());
140  tower->add_eshower(g4hit_i->get_shower_id(), g4hit_i->get_edep());
141  }
142 
143  float towerE = 0.;
144 
145  if (Verbosity())
146  {
147  towerE = m_Towers->getTotalEdep();
148  std::cout << "towers before compression: "<< m_Towers->size() << "\t" << m_Detector << std::endl;
149  }
151  if (Verbosity())
152  {
153  std::cout << "storing towers: "<< m_Towers->size() << std::endl;
154  cout << "Energy lost by dropping towers with less than " << m_Emin
155  << " energy, lost energy: " << towerE - m_Towers->getTotalEdep() << endl;
156  m_Towers->identify();
159  for (iter = begin_end.first; iter != begin_end.second; ++iter)
160  {
161  iter->second->identify();
162  }
163  }
164 
166 }
167 
169 {
171 }
172 
174 {
175  m_Detector = d;
177 }
178 
180 {
181  PHNodeIterator iter(topNode);
182  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
183  if (!runNode)
184  {
185  std::cerr << PHWHERE << "Run Node missing, doing nothing." << std::endl;
186  throw std::runtime_error("Failed to find Run node in BwdRawTowerBuilderByHitIndex::CreateNodes");
187  }
188 
189  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
190  if (!dstNode)
191  {
192  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
193  throw std::runtime_error("Failed to find DST node in BwdRawTowerBuilderByHitIndex::CreateNodes");
194  }
195 
196  // Create the tower geometry node on the tree
198  string NodeNameTowerGeometries = "TOWERGEOM_" + m_Detector;
199 
200  PHIODataNode<PHObject> *geomNode = new PHIODataNode<PHObject>(m_Geoms, NodeNameTowerGeometries, "PHObject");
201  runNode->addNode(geomNode);
202 
203  // Find detector node (or create new one if not found)
204  PHNodeIterator dstiter(dstNode);
205  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
206  "PHCompositeNode", m_Detector));
207  if (!DetNode)
208  {
209  DetNode = new PHCompositeNode(m_Detector);
210  dstNode->addNode(DetNode);
211  }
212 
213  // Create the tower nodes on the tree
215  string NodeNameTowers;
216  if (m_SimTowerNodePrefix.empty())
217  {
218  // no prefix, consistent with older convension
219  NodeNameTowers = "TOWER_" + m_Detector;
220  }
221  else
222  {
223  NodeNameTowers = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
224  }
225 
226  PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_Towers, NodeNameTowers, "PHObject");
227  DetNode->addNode(towerNode);
228 
229  return;
230 }
231 
233 {
234  /* Stream to read table from file */
235  ifstream istream_mapping;
236 
237  /* Open the datafile, if it won't open return an error */
238  if (!istream_mapping.is_open())
239  {
240  istream_mapping.open(m_MappingTowerFile);
241  if (!istream_mapping)
242  {
243  cerr << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << m_MappingTowerFile << endl;
244  exit(1);
245  }
246  }
247 
248  string line_mapping;
249 
250  while (getline(istream_mapping, line_mapping))
251  {
252  /* Skip lines starting with / including a '#' */
253  if (line_mapping.find("#") != string::npos)
254  {
255  if (Verbosity() > 0)
256  {
257  cout << "BwdRawTowerBuilderByHitIndex: SKIPPING line in mapping file: " << line_mapping << endl;
258  }
259  continue;
260  }
261 
262 // istringstream iss(line_mapping);
263 
264  /* If line starts with keyword Tower, add to tower positions */
265 // if (line_mapping.find("Tower ") != string::npos)
266 // {
267 // unsigned idx_j, idx_k, idx_l;
268 // double pos_x, pos_y, pos_z;
269 // double size_x, size_y, size_z;
270 // double rot_x, rot_y, rot_z;
271 // double type;
272 // string dummys;
273 
274  /* read string- break if error */
275 // 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))
276 // {
277 // cerr << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << endl;
278 // exit(1);
279 // }
280 
281  std::istringstream iss(line_mapping);
282 
283 // /* If line starts with keyword Tower, add to tower positions */
284 // if (line_mapping.find("Tower ") != std::string::npos)
285 // {
286  unsigned idx_j, idx_k, idx_l;
287  G4double pos_x, pos_y, pos_z, pos_z1;
288  G4double size_x, size_y, size_z;
289 // size_x = 2 *cm;
290 // size_y = 2 * cm;
291 // size_z = 20 * cm;
292  G4double rot_x, rot_y, rot_z;
293  G4double dummy;
294  std::string dummys;
295  if (line_mapping.find("Bwd ") != std::string::npos)
296  {
297  if (!(iss >> dummys>> pos_x >> pos_y >> pos_z >> pos_z1 ))
298  {
299  std::cout << "ERROR in BwdRawTowerBuilder: Failed to read global position from mapping file " << m_MappingTowerFile << std::endl;
300  gSystem->Exit(1);
301  }
302  m_GlobalPlaceInX=pos_x * cm;
303  m_GlobalPlaceInY=pos_y * cm;
304  m_GlobalPlaceInZ=(pos_z + pos_z1) * cm;
305  }
306  else if (line_mapping.find("Tower ") != std::string::npos){
307  /* read string- break if error */
308  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))
309 // if (!(iss >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> rot_x >> rot_y >> rot_z >> dummy))
310  {
311  std::cout << "ERROR in BwdRawTowerBuilder: Failed to read line in mapping file " << m_MappingTowerFile << std::endl;
312  gSystem->Exit(1);
313  }
314 
315  /* Construct unique Tower ID */
316  unsigned int temp_id = RawTowerDefs::encode_towerid(m_CaloId, idx_j, idx_k);
317  size_x = size_x * cm;
318  size_y = size_y * cm;
319  size_z = size_z * cm;
320  pos_x = pos_x * cm;
321  pos_y = pos_y * cm;
322  pos_z = 0;//48 * cm; //!!! DO NOT move tower to B0 position wrt the B0 magnet
323  /* Create tower geometry object */
324  RawTowerGeom *temp_geo = new RawTowerGeomv3(temp_id);
325  temp_geo->set_center_x(pos_x);
326  temp_geo->set_center_y(pos_y);
327  temp_geo->set_center_z(pos_z);
328  temp_geo->set_size_x(size_x);
329  temp_geo->set_size_y(size_y);
330  temp_geo->set_size_z(size_z);
331 // temp_geo->set_tower_type((int) type);
332 
333  /* Insert this tower into position map */
334  m_Geoms->add_tower_geometry(temp_geo);
335  }
336  else std::cout <<"ERROR in BwdRawTowerBuilder: 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 string / value. */
342 // string parname;
343 // double parval;
344 
345  /* read string- break if error */
346 // if (!(iss >> parname >> parval))
347 // {
348 // cerr << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << 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  cout << "* Local tower x y z : " << x_temp << " " << y_temp << " " << z_temp << endl;
388  cout << "* Tower size x y z : " << it->second->get_size_x() << " " << it->second->get_size_y() << " " << it->second->get_size_z() << endl;
389  cout << "* Globl tower x y z : " << x_temp_rt << " " << y_temp_rt << " " << z_temp_rt << endl;
390  }
391  }
392 
393  if (Verbosity())
394  {
395  cout << "size tower geom container:" << m_Geoms->size() << "\t" << m_Detector << endl;
396  }
397  return true;
398 }