ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TowerJetInput.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TowerJetInput.cc
1 #include "TowerJetInput.h"
2 
3 #include "Jet.h"
4 #include "Jetv1.h"
5 
6 #include <calobase/RawTower.h>
7 #include <calobase/RawTowerContainer.h>
8 #include <calobase/RawTowerGeom.h>
9 #include <calobase/RawTowerGeomContainer.h>
10 
11 #include <g4vertex/GlobalVertex.h>
13 
14 #include <phool/getClass.h>
15 
16 #include <cassert>
17 #include <cmath> // for asinh, atan2, cos, cosh
18 #include <iostream>
19 #include <map> // for _Rb_tree_const_iterator
20 #include <utility> // for pair
21 #include <vector>
22 
23 using namespace std;
24 
26  : _input(input)
27 {
28 }
29 
30 void TowerJetInput::identify(std::ostream &os)
31 {
32  os << " TowerJetInput: ";
33  if (_input == Jet::CEMC_TOWER)
34  os << "TOWER_CEMC to Jet::CEMC_TOWER";
35  if (_input == Jet::EEMC_TOWER)
36  os << "TOWER_EEMC to Jet::EEMC_TOWER";
37  else if (_input == Jet::HCALIN_TOWER)
38  os << "TOWER_HCALIN to Jet::HCALIN_TOWER";
39  else if (_input == Jet::HCALOUT_TOWER)
40  os << "TOWER_HCALOUT to Jet::HCALOUT_TOWER";
41  else if (_input == Jet::FEMC_TOWER)
42  os << "TOWER_FEMC to Jet::FEMC_TOWER";
43  else if (_input == Jet::FHCAL_TOWER)
44  os << "TOWER_FHCAL to Jet::FHCAL_TOWER";
45  os << endl;
46 }
47 
48 std::vector<Jet *> TowerJetInput::get_input(PHCompositeNode *topNode)
49 {
50  if (Verbosity() > 0) cout << "TowerJetInput::process_event -- entered" << endl;
51 
52  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
53  if (!vertexmap)
54  {
55  cout << "TowerJetInput::get_input - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
56  assert(vertexmap); // force quit
57 
58  return std::vector<Jet *>();
59  }
60 
61  if (vertexmap->empty())
62  {
63  cout << "TowerJetInput::get_input - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_bbc or tracking reco flags in the main macro in order to reconstruct the global vertex." << endl;
64  return std::vector<Jet *>();
65  }
66 
67  RawTowerContainer *towers = nullptr;
68  RawTowerGeomContainer *geom = nullptr;
69  if (_input == Jet::CEMC_TOWER)
70  {
71  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
72  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
73  if (!towers || !geom)
74  {
75  return std::vector<Jet *>();
76  }
77  }
78  else if (_input == Jet::EEMC_TOWER)
79  {
80  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_EEMC");
81  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_EEMC");
82  if (!towers || !geom)
83  {
84  return std::vector<Jet *>();
85  }
86  }
87  else if (_input == Jet::HCALIN_TOWER)
88  {
89  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
90  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
91  if (!towers || !geom)
92  {
93  return std::vector<Jet *>();
94  }
95  }
96  else if (_input == Jet::HCALOUT_TOWER)
97  {
98  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
99  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
100  if (!towers || !geom)
101  {
102  return std::vector<Jet *>();
103  }
104  }
105  else if (_input == Jet::FEMC_TOWER)
106  {
107  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_FEMC");
108  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_FEMC");
109  if (!towers || !geom)
110  {
111  return std::vector<Jet *>();
112  }
113  }
114  else if (_input == Jet::FHCAL_TOWER)
115  {
116  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_FHCAL");
117  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_FHCAL");
118  if (!towers || !geom)
119  {
120  return std::vector<Jet *>();
121  }
122  }
123  else if (_input == Jet::CEMC_TOWER_RETOWER)
124  {
125  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
126  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
127  if (!towers || !geom)
128  {
129  return std::vector<Jet *>();
130  }
131  }
132  else if (_input == Jet::CEMC_TOWER_SUB1)
133  {
134  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
135  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
136  if (!towers || !geom)
137  {
138  return std::vector<Jet *>();
139  }
140  }
141  else if (_input == Jet::HCALIN_TOWER_SUB1)
142  {
143  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
144  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
145  if (!towers || !geom)
146  {
147  return std::vector<Jet *>();
148  }
149  }
150  else if (_input == Jet::HCALOUT_TOWER_SUB1)
151  {
152  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
153  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
154  if (!towers || !geom)
155  {
156  return std::vector<Jet *>();
157  }
158  }
159  else if (_input == Jet::CEMC_TOWER_SUB1CS)
160  {
161  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1CS");
162  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
163  if (!towers || !geom)
164  {
165  return std::vector<Jet *>();
166  }
167  }
168  else if (_input == Jet::HCALIN_TOWER_SUB1CS)
169  {
170  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1CS");
171  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
172  if (!towers || !geom)
173  {
174  return std::vector<Jet *>();
175  }
176  }
177  else if (_input == Jet::HCALOUT_TOWER_SUB1CS)
178  {
179  towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1CS");
180  geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
181  if (!towers || !geom)
182  {
183  return std::vector<Jet *>();
184  }
185  }
186  else
187  {
188  return std::vector<Jet *>();
189  }
190 
191  // first grab the event vertex or bail
192  GlobalVertex *vtx = vertexmap->begin()->second;
193  float vtxz = NAN;
194  if (vtx)
195  vtxz = vtx->get_z();
196  else
197  return std::vector<Jet *>();
198 
199  if (isnan(vtxz))
200  {
201  static bool once = true;
202  if (once)
203  {
204  once = false;
205 
206  cout << "TowerJetInput::get_input - WARNING - vertex is NAN. Drop all tower inputs (further NAN-vertex warning will be suppressed)." << endl;
207  }
208 
209  return std::vector<Jet *>();
210  }
211 
212  std::vector<Jet *> pseudojets;
213  RawTowerContainer::ConstRange begin_end = towers->getTowers();
215  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
216  {
217  RawTower *tower = rtiter->second;
218 
219  RawTowerGeom *tower_geom =
220  geom->get_tower_geometry(tower->get_key());
221  assert(tower_geom);
222 
223  double r = tower_geom->get_center_radius();
224  double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
225  double z0 = tower_geom->get_center_z();
226 
227  double z = z0 - vtxz;
228 
229  double eta = asinh(z / r); // eta after shift from vertex
230 
231  double pt = tower->get_energy() / cosh(eta);
232  double px = pt * cos(phi);
233  double py = pt * sin(phi);
234  double pz = pt * sinh(eta);
235 
236  Jet *jet = new Jetv1();
237  jet->set_px(px);
238  jet->set_py(py);
239  jet->set_pz(pz);
240  jet->set_e(tower->get_energy());
241  jet->insert_comp(_input, tower->get_id());
242 
243  pseudojets.push_back(jet);
244  }
245 
246  if (Verbosity() > 0) cout << "TowerJetInput::process_event -- exited" << endl;
247 
248  return pseudojets;
249 }