ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerBuilderDRCALO.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerBuilderDRCALO.cc
2 
3 #include <Geant4/G4SystemOfUnits.hh>
4 #include <Geant4/G4Types.hh> // for G4double, G4int
5 
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerv1.h>
8 #include <calobase/RawTowerv2.h>
9 
10 #include <calobase/RawTower.h> // for RawTower
11 #include <calobase/RawTowerDefs.h> // for convert_name_to_caloid
12 #include <calobase/RawTowerGeom.h> // for RawTowerGeom
13 #include <calobase/RawTowerGeomContainer.h> // for RawTowerGeomContainer
14 #include <calobase/RawTowerGeomContainerv1.h>
15 #include <calobase/RawTowerGeomv3.h>
16 
17 #include <g4main/PHG4Hit.h>
19 
21 #include <fun4all/SubsysReco.h> // for SubsysReco
22 
23 #include <phool/PHCompositeNode.h>
24 #include <phool/PHIODataNode.h>
25 #include <phool/PHNode.h> // for PHNode
26 #include <phool/PHNodeIterator.h>
27 #include <phool/PHObject.h> // for PHObject
28 #include <phool/getClass.h>
29 #include <phool/phool.h> // for PHWHERE
30 
31 #include <TRotation.h>
32 #include <TVector3.h>
33 
34 #include <cstdlib> // for exit
35 #include <exception> // for exception
36 #include <fstream>
37 #include <iostream>
38 #include <map>
39 #include <sstream>
40 #include <stdexcept>
41 #include <utility> // for pair, make_pair
42 
43 using namespace std;
44 
46  : SubsysReco(name)
47  , m_Towers(nullptr)
48  , m_Geoms(nullptr)
49  , m_Detector("NONE")
50  , m_MappingTowerFile("default.txt")
51  , m_CaloId(RawTowerDefs::NONE)
52  , m_GlobalPlaceInX(0)
53  , m_GlobalPlaceInY(0)
54  , m_GlobalPlaceInZ(0)
55  , m_Emin(1e-6)
56 {
57 }
58 
60 {
61  PHNodeIterator iter(topNode);
62 
63  // Looking for the DST node
64  PHCompositeNode *dstNode;
65  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
66  if (!dstNode)
67  {
68  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
69  exit(1);
70  }
71 
72  try
73  {
74  CreateNodes(topNode);
75  }
76  catch (std::exception &e)
77  {
78  std::cout << e.what() << std::endl;
79  //exit(1);
80  }
81 
82  try
83  {
85  }
86  catch (std::exception &e)
87  {
88  std::cout << e.what() << std::endl;
89  //exit(1);
90  }
91 
93 }
94 
96 {
97  // get hits
98  string NodeNameHits = "G4HIT_" + m_Detector;
99  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, NodeNameHits);
100  if (!g4hit)
101  {
102  cout << "Could not locate g4 hit node " << NodeNameHits << endl;
103  exit(1);
104  }
105 
106  // loop over all hits in the event
108  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
109 
110  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
111  {
112  PHG4Hit *g4hit_i = hiter->second;
113 
114  // Don't include hits with zero energy
115  if (g4hit_i->get_edep() <= 0 && g4hit_i->get_edep() != -1) continue;
116  // cout << g4hit_i->get_index_j() << "\t" << g4hit_i->get_index_k() << "\t" << endl;
117  // encode CaloTowerID from j, k index of tower / hit and calorimeter ID
119  g4hit_i->get_index_j(),
120  g4hit_i->get_index_k());
121  // add the energy to the corresponding tower
122  RawTowerv2 *tower = dynamic_cast<RawTowerv2 *>(m_Towers->getTower(calotowerid));
123  if (!tower)
124  {
125  tower = new RawTowerv2(calotowerid);
126  tower->set_energy(0);
127  tower->set_scint_gammas(0.);
128  tower->set_cerenkov_gammas(0.);
129  m_Towers->AddTower(tower->get_id(), tower);
130 // cout << "intializing tower" << endl;
131 // tower->identify();
132  }
133  else
134  {
135  // check version consistency
136  if (dynamic_cast<RawTowerv2 *>(tower) == nullptr)
137  {
138  cout << __PRETTY_FUNCTION__ << " : Fatal Error! "
139  << "Expect RawTowerv2, but found this tower:";
140  tower->identify();
141  exit(1);
142  }
143  }
144 
145  if (Verbosity() > 3)
146  cout << g4hit_i->get_property_float(PHG4Hit::PROPERTY::scint_gammas) << "\t" << g4hit_i->get_property_float(PHG4Hit::PROPERTY::cerenkov_gammas) << endl;
147  // cout << tower->get_scint_gammas() << "\t" << g4hit_i->get_property_float(PHG4Hit::PROPERTY::scint_gammas) << "\t" << tower->get_cerenkov_gammas() << "\t" << g4hit_i->get_property_float(PHG4Hit::PROPERTY::cerenkov_gammas) << endl;
148  tower->set_scint_gammas(tower->get_scint_gammas() + g4hit_i->get_property_float(PHG4Hit::PROPERTY::scint_gammas));
149  tower->set_cerenkov_gammas(tower->get_cerenkov_gammas() + g4hit_i->get_property_float(PHG4Hit::PROPERTY::cerenkov_gammas));
150 
151  tower->add_ecell((g4hit_i->get_index_j() << 16) + g4hit_i->get_index_k(), g4hit_i->get_light_yield());
152  // cout << (g4hit_i->get_index_j() << 16) + g4hit_i->get_index_k() << "\t" << g4hit_i->get_light_yield() << endl;
153  tower->set_energy(tower->get_energy() + g4hit_i->get_light_yield());
154  tower->add_eshower(g4hit_i->get_shower_id(), g4hit_i->get_edep());
155  // tower->identify();
156  }
157 
158  float towerE = 0.;
159 
160  if (Verbosity())
161  {
162  towerE = m_Towers->getTotalEdep();
163  }
164 
166  if (Verbosity())
167  {
168  cout << "Energy lost by dropping towers with less than " << m_Emin
169  << " energy, lost energy: " << towerE - m_Towers->getTotalEdep() << endl;
170  m_Towers->identify();
173  for (iter = begin_end.first; iter != begin_end.second; ++iter)
174  {
175  iter->second->identify();
176  }
177  }
178 
180 }
181 
183 {
185 }
186 
187 void RawTowerBuilderDRCALO::Detector(const std::string &d)
188 {
189  m_Detector = d;
191 }
192 
194 {
195  PHNodeIterator iter(topNode);
196  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
197  if (!runNode)
198  {
199  std::cerr << PHWHERE << "Run Node missing, doing nothing." << std::endl;
200  throw std::runtime_error("Failed to find Run node in RawTowerBuilderDRCALO::CreateNodes");
201  }
202 
203  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
204  if (!dstNode)
205  {
206  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
207  throw std::runtime_error("Failed to find DST node in RawTowerBuilderDRCALO::CreateNodes");
208  }
209 
210  // Create the tower geometry node on the tree
212  string NodeNameTowerGeometries = "TOWERGEOM_" + m_Detector;
213 
214  PHIODataNode<PHObject> *geomNode = new PHIODataNode<PHObject>(m_Geoms, NodeNameTowerGeometries, "PHObject");
215  runNode->addNode(geomNode);
216 
217  // Find detector node (or create new one if not found)
218  PHNodeIterator dstiter(dstNode);
219  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
220  "PHCompositeNode", m_Detector));
221  if (!DetNode)
222  {
223  DetNode = new PHCompositeNode(m_Detector);
224  dstNode->addNode(DetNode);
225  }
226 
227  // Create the tower nodes on the tree
229  string NodeNameTowers;
230  if (m_SimTowerNodePrefix.empty())
231  {
232  // no prefix, consistent with older convension
233  NodeNameTowers = "TOWER_" + m_Detector;
234  }
235  else
236  {
237  NodeNameTowers = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
238  }
239 
240  PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_Towers, NodeNameTowers, "PHObject");
241  DetNode->addNode(towerNode);
242 
243  return;
244 }
245 // /*
247 {
248  // // Stream to read table from file
249  ifstream istream_mapping;
250 
251  // // Open the datafile, if it won't open return an error
252  if (!istream_mapping.is_open())
253  {
254  istream_mapping.open(m_MappingTowerFile);
255  if (!istream_mapping)
256  {
257  cerr << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << m_MappingTowerFile << endl;
258  exit(1);
259  }
260  }
261 
262  string line_mapping;
263 
264  while (getline(istream_mapping, line_mapping))
265  {
266  // // Skip lines starting with / including a '#'
267  if (line_mapping.find("#") != string::npos)
268  {
269  if (Verbosity() > 0)
270  {
271  cout << "RawTowerBuilderDRCALO: SKIPPING line in mapping file: " << line_mapping << endl;
272  }
273  continue;
274  }
275 
276  istringstream iss(line_mapping);
277 
278  // // If line starts with keyword Tower, add to tower positions
279  if (line_mapping.find("Tower ") != string::npos)
280  {
281  }
282 
283  // // If line does not start with keyword Tower, read as parameter
284  else
285  {
286  // // If this line is not a comment and not a tower, save parameter as string / value.
287  string parname;
288  double parval;
289 
290  // // read string- break if error
291  if (!(iss >> parname >> parval))
292  {
293  cerr << "ERROR in RawTowerBuilderDRCALO: Failed to read line in mapping file " << m_MappingTowerFile << endl;
294  exit(1);
295  }
296 
297  m_GlobalParameterMap.insert(make_pair(parname, parval));
298  }
299  }
300 
301  // // Update member variables for global parameters based on parsed parameter file
302  std::map<string, double>::iterator parit;
303 
304  parit = m_GlobalParameterMap.find("Gx0");
305  if (parit != m_GlobalParameterMap.end())
306  m_GlobalPlaceInX = parit->second * cm;
307 
308  parit = m_GlobalParameterMap.find("Gy0");
309  if (parit != m_GlobalParameterMap.end())
310  m_GlobalPlaceInY = parit->second * cm;
311 
312  parit = m_GlobalParameterMap.find("Gz0");
313  if (parit != m_GlobalParameterMap.end())
314  m_GlobalPlaceInZ = parit->second * cm;
315 
316  // unsigned idx_j, idx_k, idx_l;
317  G4double pos_x, pos_y, pos_z;
318  // double size_x, size_y, size_z;
319  G4double size_z = 10;
320  // double rot_x, rot_y, rot_z;
321  double type = 0;
322  string dummys;
323 
324  parit = m_GlobalParameterMap.find("Gdz");
325  if (parit != m_GlobalParameterMap.end())
326  size_z = parit->second * cm;
327 
328  pos_z = m_GlobalPlaceInZ;
329  G4double twredge = 1.0;//1.2
330  G4double twrreadout = 1.0;//1.2
331  G4double twrsize = 1.0;//1.2
332  G4double drsize = 220;
333  parit = m_GlobalParameterMap.find("Gtower_dx");
334  if (parit != m_GlobalParameterMap.end()){
335  twrsize = parit->second * cm;
336  }
337 
338  parit = m_GlobalParameterMap.find("Gtower_readout");
339  if (parit != m_GlobalParameterMap.end())
340  {
341  twrreadout = parit->second * cm;
342  }
343  parit = m_GlobalParameterMap.find("Gr1_outer");
344  if (parit != m_GlobalParameterMap.end()){
345  drsize = parit->second * cm;
346  }
347  int maxsubtow = (int) ( (twrsize) / (twrreadout));
348  twredge = (twrsize - (maxsubtow * twrreadout))/maxsubtow;
349  if(twrsize!=twrreadout){
350  twrsize = twrreadout + twredge;
351  }
352 
353  float scaling = 1.;
354  for (int idxj = 0; idxj < drsize * 2 / twrsize; idxj++)
355  {
356  pos_x = (idxj * twrsize - drsize); //TODO DRCALO TOWER SIZE
357  // cout << pos_x << endl;
358  for (int idxk = 0; idxk < drsize * 2 / twrsize; idxk++)
359  {
360  pos_y = (idxk * twrsize - drsize); //TODO DRCALO TOWER SIZE
361  // // Construct unique Tower ID
362  unsigned int temp_id = RawTowerDefs::encode_towerid(m_CaloId, idxj, idxk);
363 
364  // // Create tower geometry object
365  RawTowerGeom *temp_geo = new RawTowerGeomv3(temp_id);
366  temp_geo->set_center_x(scaling * pos_x);
367  temp_geo->set_center_y(scaling * pos_y);
368  temp_geo->set_center_z(pos_z);
369  temp_geo->set_size_x(scaling * twrsize); //TODO DRCALO TOWER SIZE
370  temp_geo->set_size_y(scaling * twrsize); //TODO DRCALO TOWER SIZE
371  temp_geo->set_size_z(size_z);
372  temp_geo->set_tower_type((int) type);
373 
374  // // Insert this tower into position map
375  m_Geoms->add_tower_geometry(temp_geo);
376  }
377  }
378 
379  return true;
380 }
381 // */