Analysis Software
Documentation for sPHENIX simulation 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/TowerInfo.h>
16 #include <calobase/TowerInfoContainerv1.h>
17 #include <calobase/RawTowerGeom.h>
18 #include <calobase/RawTowerGeomContainer.h>
19 
22 
23 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY_MORE
24 #include <fun4all/SubsysReco.h>
25 
26 #include <phool/getClass.h>
27 
28 #include <CLHEP/Vector/ThreeVector.h>
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 
118 /*const*/ CLHEP::Hep3Vector ClusterIso::getVertex()
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  TowerInfoContainer *towersEM3old = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERINFO_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() << " TOWERINFO_CALIB_CEMC_RETOWER_SUB1 towers" << '\n';
151 
152  //get InnerHCal towers
153  TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
154  if (Verbosity() >= VERBOSITY_MORE) std::cout << Name() << "::ClusterIso::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN_SUB1 towers" << '\n';
155 
156  //get outerHCal towers
157  TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
158  if (Verbosity() >= VERBOSITY_MORE) std::cout << Name() << "::ClusterIso::process_event: " << towersOH3->size() << " TOWERINFO_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  unsigned int ntowers = towersEM3old->size();
206  for (unsigned int channel = 0; channel < ntowers;channel++)
207  {
208  TowerInfo *tower = towersEM3old->get_tower_at_channel(channel);
209  unsigned int towerkey = towersEM3old->encode_key(channel);
210  int ieta = towersEM3old->getTowerEtaBin(towerkey);
211  int iphi = towersEM3old->getTowerPhiBin(towerkey);
212 
214 
215  RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
216  double this_phi = tower_geom->get_phi();
217  double this_eta = tower_geom->get_eta();
218  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
219  {
220  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
221  }
222  }
223  }
224 
225  //calculate Inner HCal tower contribution to isolation energy
226  {
227  unsigned int ntowers = towersIH3->size();
228  for (unsigned int channel = 0; channel < ntowers;channel++)
229  {
231  unsigned int towerkey =towersIH3->encode_key(channel);
232  int ieta = towersIH3->getTowerEtaBin(towerkey);
233  int iphi = towersIH3->getTowerPhiBin(towerkey);
235  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(key);
236  double this_phi = tower_geom->get_phi();
237  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
238  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
239  {
240  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
241  }
242  }
243  }
244 
245  //calculate Outer HCal tower contribution to isolation energy
246  {
247  unsigned int ntowers = towersOH3->size();
248  for (unsigned int channel = 0; channel < ntowers;channel++)
249  {
251  unsigned int towerkey =towersOH3->encode_key(channel);
252  int ieta = towersOH3->getTowerEtaBin(towerkey);
253  int iphi = towersOH3->getTowerPhiBin(towerkey);
255  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(key);
256  double this_phi = tower_geom->get_phi();
257  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
258  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
259  {
260  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
261  }
262  }
263  }
264 
265  isoEt -= et; //Subtract cluster eT from isoET
267  {
268  std::cout << Name() << "::ClusterIso iso_et for ";
269  cluster->identify();
270  std::cout << "=" << isoEt << '\n';
271  }
272  cluster->set_et_iso(isoEt, (int) 10 * m_coneSize, 1, 1);
273  }
274  }
275  }
276  }
277  if (m_do_unsubtracted)
278  {
282  if (Verbosity() >= VERBOSITY_EVEN_MORE) std::cout << Name() << "::ClusterIso starting unsubtracted calculation" << '\n';
283  {
284  //get EMCal towers
285  TowerInfoContainer *towersEM3old = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERINFO_CALIB_CEMC");
286  if (Verbosity() >= VERBOSITY_MORE) std::cout << "ClusterIso::process_event: " << towersEM3old->size() << " TOWERINFO_CALIB_CEMC towers" << '\n';
287 
288  //get InnerHCal towers
289  TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERINFO_CALIB_HCALIN");
290  if (Verbosity() >= VERBOSITY_MORE) std::cout << "ClusterIso::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN towers" << '\n';
291 
292  //get outerHCal towers
293  TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERINFO_CALIB_HCALOUT");
294  if (Verbosity() >= VERBOSITY_MORE) std::cout << "ClusterIso::process_event: " << towersOH3->size() << " TOWERINFO_CALIB_HCALOUT towers" << std::endl;
295 
296  //get geometry of calorimeter towers
297  RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
298  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
299  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
300 
301  {
302  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
303  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
305  if (Verbosity() >= VERBOSITY_SOME) std::cout << " ClusterIso sees " << clusters->size() << " clusters " << '\n';
306 
307  //vertexmap is used to get correct collision vertex
308  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
309  m_vx = m_vy = m_vz = 0;
310  if (vertexmap && !vertexmap->empty())
311  {
312  GlobalVertex *vertex = (vertexmap->begin()->second);
313  m_vx = vertex->get_x();
314  m_vy = vertex->get_y();
315  m_vz = vertex->get_z();
316  if (Verbosity() >= VERBOSITY_SOME) std::cout << Name() << "ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
317  }
318 
319  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
320  {
321  RawCluster *cluster = rtiter->second;
322 
323  CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
324  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
325  double cluster_energy = E_vec_cluster.mag();
326  double cluster_eta = E_vec_cluster.pseudoRapidity();
327  double cluster_phi = E_vec_cluster.phi();
328  double et = cluster_energy / cosh(cluster_eta);
329  double isoEt = 0;
330  if (Verbosity() >= VERBOSITY_MAX)
331  {
332  std::cout << Name() << "::ClusterIso processing";
333  cluster->identify();
334  std::cout << '\n';
335  }
336  if (et < m_eTCut)
337  {
338  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t does not pass eT cut" << '\n';
339  continue;
340  } //skip if cluster is below eT cut
341 
342  //calculate EMCal tower contribution to isolation energy
343  {
344  unsigned int ntowers = towersEM3old->size();
345  for (unsigned int channel = 0; channel < ntowers;channel++)
346  {
347  TowerInfo *tower = towersEM3old->get_tower_at_channel(channel);
348  unsigned int towerkey = towersEM3old->encode_key(channel);
349  int ieta = towersEM3old->getTowerEtaBin(towerkey);
350  int iphi = towersEM3old->getTowerPhiBin(towerkey);
352  RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
353  double this_phi = tower_geom->get_phi();
354  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
355  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
356  {
357  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
358  }
359  }
360  }
361  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t after EMCal isoEt:" << isoEt << '\n';
362  //calculate Inner HCal tower contribution to isolation energy
363  {
364  unsigned int ntowers = towersIH3->size();
365  for (unsigned int channel = 0; channel < ntowers;channel++)
366  {
368  unsigned int towerkey = towersIH3->encode_key(channel);
369  int ieta = towersIH3->getTowerEtaBin(towerkey);
370  int iphi = towersIH3->getTowerPhiBin(towerkey);
372  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(key);
373  double this_phi = tower_geom->get_phi();
374  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
375  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
376  {
377  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
378  }
379  }
380  }
381  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t after innerHCal isoEt:" << isoEt << '\n';
382  //calculate Outer HCal tower contribution to isolation energy
383  {
384  unsigned int ntowers = towersOH3->size();
385  for (unsigned int channel = 0; channel < ntowers;channel++)
386  {
388  unsigned int towerkey = towersOH3->encode_key(channel);
389  int ieta = towersOH3->getTowerEtaBin(towerkey);
390  int iphi = towersOH3->getTowerPhiBin(towerkey);
392  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(key);
393  double this_phi = tower_geom->get_phi();
394  double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
395  if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
396  {
397  isoEt += tower->get_energy() / cosh(this_eta); //if tower is in cone, add energy
398  }
399  }
400  }
401  if (Verbosity() >= VERBOSITY_MAX) std::cout << "\t after outerHCal isoEt:" << isoEt << '\n';
402  isoEt -= et; //Subtract cluster eT from isoET
404  {
405  std::cout << Name() << "::ClusterIso iso_et for ";
406  cluster->identify();
407  std::cout << "=" << isoEt << '\n';
408  }
409  cluster->set_et_iso(isoEt, (int) 10 * m_coneSize, 0, 1);
410  }
411  }
412  }
413  }
414  return 0;
415 }
416 
418 {
419  return 0;
420 }