ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ClusterIso.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ClusterIso.cc
1 
10 #include "ClusterIso.h"
11 
12 #include <calobase/RawCluster.h>
13 #include <calobase/RawClusterContainer.h>
14 #include <calobase/RawClusterUtility.h>
15 #include <calobase/RawTower.h>
16 #include <calobase/RawTowerContainer.h>
17 #include <calobase/RawTowerGeom.h>
18 #include <calobase/RawTowerGeomContainer.h>
19 
20 #include <g4vertex/GlobalVertex.h>
22 
23 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY_MORE
24 #include <fun4all/SubsysReco.h>
25 
26 #include <phool/getClass.h>
27 
29 
30 #include <iostream>
31 #include <map>
32 #include <utility>
33 
40 double ClusterIso::getTowerEta(RawTowerGeom *tower_geom, double vx, double vy, double vz)
41 {
42  float r;
43  if (vx == 0 && vy == 0 && vz == 0)
44  {
45  r = tower_geom->get_eta();
46  }
47  else
48  {
49  double radius = sqrt((tower_geom->get_center_x() - vx) * (tower_geom->get_center_x() - vx) + (tower_geom->get_center_y() - vy) * (tower_geom->get_center_y() - vy));
50  double theta = atan2(radius, tower_geom->get_center_z() - vz);
51  r = -log(tan(theta / 2.));
52  }
53  return r;
54 }
55 
60 ClusterIso::ClusterIso(const std::string &kname, float eTCut = 0.0, int coneSize = 3, bool do_subtracted = 1, bool do_unsubtracted = 1)
61  : SubsysReco(kname)
62  , m_do_subtracted(do_subtracted)
63  , m_do_unsubtracted(do_unsubtracted)
64 {
65  if (Verbosity() >= VERBOSITY_SOME) std::cout << Name() << "::ClusterIso constructed" << '\n';
66  if (coneSize == 0 && Verbosity() >= VERBOSITY_QUIET) std::cout << "WARNING in " << Name() << "ClusterIso:: cone size is zero" << '\n';
67  m_vx = m_vy = m_vz = 0;
68  setConeSize(coneSize);
69  seteTCut(eTCut);
71  {
72  std::cout << Name() << "::ClusterIso::m_coneSize is:" << m_coneSize << '\n';
73  std::cout << Name() << "::ClusterIso::m_eTCut is:" << m_eTCut << '\n';
74  }
75  if (!do_subtracted && !do_unsubtracted && Verbosity() >= VERBOSITY_QUIET) std::cout << "WARNING in " << Name() << "ClusterIso:: all processes turned off doing nothing" << '\n';
76 }
77 
79 {
80  return 0;
81 }
82 
86 void ClusterIso::seteTCut(float eTCut)
87 {
88  this->m_eTCut = eTCut;
89 }
90 
94 void ClusterIso::setConeSize(int coneSize)
95 {
96  this->m_coneSize = coneSize / 10.0;
97 }
98 
102 /*const*/ float ClusterIso::geteTCut()
103 {
104  return m_eTCut;
105 }
106 
110 /*const*/ int ClusterIso::getConeSize()
111 {
112  return (int) m_coneSize * 10;
113 }
114 
119 {
120  return CLHEP::Hep3Vector(m_vx, m_vy, m_vz);
121 }
122 
130 {
131  if (Verbosity() >= VERBOSITY_MORE) std::cout << Name() << "::ClusterIso::process_event" << '\n';
139  if (m_do_subtracted)
140  {
141  {
142  if (Verbosity() >= VERBOSITY_EVEN_MORE) std::cout << Name() << "::ClusterIso starting subtracted calculation" << '\n';
143  //get EMCal towers
144  RawTowerContainer *towersEM3old = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
145  if (towersEM3old == nullptr)
146  {
147  m_do_subtracted = false;
148  if (Verbosity() >= VERBOSITY_SOME) std::cout << "In " << Name() << "::ClusterIso WARNING substracted towers do not exist subtracted isolation cannot be preformed \n";
149  }
150  if (Verbosity() >= VERBOSITY_MORE) std::cout << Name() << "::ClusterIso::process_event: " << towersEM3old->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1 towers" << '\n';
151 
152  //get InnerHCal towers
153  RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
154  if (Verbosity() >= VERBOSITY_MORE) std::cout << Name() << "::ClusterIso::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN_SUB1 towers" << '\n';
155 
156  //get outerHCal towers
157  RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
158  if (Verbosity() >= VERBOSITY_MORE) std::cout << Name() << "::ClusterIso::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT_SUB1 towers" << std::endl;
159 
160  //get geometry of calorimeter towers
161  RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
162  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
163  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
164 
165  {
166  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
167  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
169  if (Verbosity() >= VERBOSITY_SOME) std::cout << Name() << "::ClusterIso sees " << clusters->size() << " clusters " << '\n';
170 
171  //vertexmap is used to get correct collision vertex
172  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
173  m_vx = m_vy = m_vz = 0;
174  if (vertexmap && !vertexmap->empty())
175  {
176  GlobalVertex *vertex = (vertexmap->begin()->second);
177  m_vx = vertex->get_x();
178  m_vy = vertex->get_y();
179  m_vz = vertex->get_z();
180  if (Verbosity() >= VERBOSITY_SOME)
181  {
182  std::cout << Name() << "::ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
183  }
184  }
185 
186  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
187  {
188  RawCluster *cluster = rtiter->second;
189 
190  CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
191  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
192  double cluster_energy = E_vec_cluster.mag();
193  double cluster_eta = E_vec_cluster.pseudoRapidity();
194  double cluster_phi = E_vec_cluster.phi();
195  double et = cluster_energy / cosh(cluster_eta);
196  double isoEt = 0;
197 
198  if (et < m_eTCut)
199  {
200  continue;
201  } //skip if cluster is under eT cut
202 
203  //calculate EMCal tower contribution to isolation energy
204  {
205  RawTowerContainer::ConstRange begin_end = towersEM3old->getTowers();
206  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
207  {
208  RawTower *tower = rtiter->second;
209  RawTowerGeom *tower_geom = geomEM->get_tower_geometry(tower->get_key());
210  double this_phi = tower_geom->get_phi();
211  double this_eta = tower_geom->get_eta();
212  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
213  {
214  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
215  }
216  }
217  }
218 
219  //calculate Inner HCal tower contribution to isolation energy
220  {
221  RawTowerContainer::ConstRange begin_end = towersIH3->getTowers();
222  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
223  {
224  RawTower *tower = rtiter->second;
225  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
226  double this_phi = tower_geom->get_phi();
227  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
228  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
229  {
230  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
231  }
232  }
233  }
234 
235  //calculate Outer HCal tower contribution to isolation energy
236  {
237  RawTowerContainer::ConstRange begin_end = towersOH3->getTowers();
238  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
239  {
240  RawTower *tower = rtiter->second;
241  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
242  double this_phi = tower_geom->get_phi();
243  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
244  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
245  {
246  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
247  }
248  }
249  }
250 
251  isoEt -= et; //Subtract cluster eT from isoET
253  {
254  std::cout << Name() << "::ClusterIso iso_et for ";
255  cluster->identify();
256  std::cout << "=" << isoEt << '\n';
257  }
258  cluster->set_et_iso(isoEt, (int) 10 * m_coneSize, 1, 1);
259  }
260  }
261  }
262  }
263  if (m_do_unsubtracted)
264  {
268  if (Verbosity() >= VERBOSITY_EVEN_MORE) std::cout << Name() << "::ClusterIso starting unsubtracted calculation" << '\n';
269  {
270  //get EMCal towers
271  RawTowerContainer *towersEM3old = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
272  if (Verbosity() >= VERBOSITY_MORE) std::cout << "ClusterIso::process_event: " << towersEM3old->size() << " TOWER_CALIB_CEMC towers" << '\n';
273 
274  //get InnerHCal towers
275  RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
276  if (Verbosity() >= VERBOSITY_MORE) std::cout << "ClusterIso::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << '\n';
277 
278  //get outerHCal towers
279  RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
280  if (Verbosity() >= VERBOSITY_MORE) std::cout << "ClusterIso::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
281 
282  //get geometry of calorimeter towers
283  RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
284  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
285  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
286 
287  {
288  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
289  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
291  if (Verbosity() >= VERBOSITY_SOME) std::cout << " ClusterIso sees " << clusters->size() << " clusters " << '\n';
292 
293  //vertexmap is used to get correct collision vertex
294  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
295  m_vx = m_vy = m_vz = 0;
296  if (vertexmap && !vertexmap->empty())
297  {
298  GlobalVertex *vertex = (vertexmap->begin()->second);
299  m_vx = vertex->get_x();
300  m_vy = vertex->get_y();
301  m_vz = vertex->get_z();
302  if (Verbosity() >= VERBOSITY_SOME) std::cout << Name() << "ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
303  }
304 
305  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
306  {
307  RawCluster *cluster = rtiter->second;
308 
309  CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
310  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
311  double cluster_energy = E_vec_cluster.mag();
312  double cluster_eta = E_vec_cluster.pseudoRapidity();
313  double cluster_phi = E_vec_cluster.phi();
314  double et = cluster_energy / cosh(cluster_eta);
315  double isoEt = 0;
316  if (Verbosity() >= VERBOSITY_MAX)
317  {
318  std::cout << Name() << "::ClusterIso processing";
319  cluster->identify();
320  std::cout << '\n';
321  }
322  if (et < m_eTCut)
323  {
324  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t does not pass eT cut" << '\n';
325  continue;
326  } //skip if cluster is below eT cut
327 
328  //calculate EMCal tower contribution to isolation energy
329  {
330  RawTowerContainer::ConstRange begin_end = towersEM3old->getTowers();
331  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
332  {
333  RawTower *tower = rtiter->second;
334  RawTowerGeom *tower_geom = geomEM->get_tower_geometry(tower->get_key());
335  double this_phi = tower_geom->get_phi();
336  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
337  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
338  {
339  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
340  }
341  }
342  }
343  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t after EMCal isoEt:" << isoEt << '\n';
344  //calculate Inner HCal tower contribution to isolation energy
345  {
346  RawTowerContainer::ConstRange begin_end = towersIH3->getTowers();
347  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
348  {
349  RawTower *tower = rtiter->second;
350  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
351  double this_phi = tower_geom->get_phi();
352  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
353  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
354  {
355  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
356  }
357  }
358  }
359  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t after innerHCal isoEt:" << isoEt << '\n';
360  //calculate Outer HCal tower contribution to isolation energy
361  {
362  RawTowerContainer::ConstRange begin_end = towersOH3->getTowers();
363  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
364  {
365  RawTower *tower = rtiter->second;
366  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
367  double this_phi = tower_geom->get_phi();
368  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
369  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
370  {
371  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
372  }
373  }
374  }
375  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t after outerHCal isoEt:" << isoEt << '\n';
376  isoEt -= et; //Subtract cluster eT from isoET
378  {
379  std::cout << Name() << "::ClusterIso iso_et for ";
380  cluster->identify();
381  std::cout << "=" << isoEt << '\n';
382  }
383  cluster->set_et_iso(isoEt, (int) 10 * m_coneSize, 0, 1);
384  }
385  }
386  }
387  }
388  return 0;
389 }
390 
392 {
393  return 0;
394 }