Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DetermineTowerBackground.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DetermineTowerBackground.cc
1 #include "DetermineTowerBackground.h"
2 
3 #include "TowerBackground.h"
4 #include "TowerBackgroundv1.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 <calobase/TowerInfo.h>
12 #include <calobase/TowerInfoContainer.h>
13 
14 #include <jetbase/Jet.h>
15 #include <jetbase/JetContainer.h>
16 
17 #include <g4main/PHG4Particle.h>
18 #include <g4main/PHG4TruthInfoContainer.h>
19 
20 #include <fun4all/Fun4AllReturnCodes.h>
21 #include <fun4all/SubsysReco.h>
22 
23 #include <phool/PHCompositeNode.h>
24 #include <phool/PHIODataNode.h>
25 #include <phool/PHNode.h>
26 #include <phool/PHNodeIterator.h>
27 #include <phool/PHObject.h>
28 #include <phool/getClass.h>
29 #include <phool/phool.h>
30 
31 #include <TLorentzVector.h>
32 
33 // standard includes
34 #include <cmath>
35 #include <cstdlib>
36 #include <iomanip>
37 #include <iostream>
38 #include <map>
39 #include <memory>
40 #include <vector>
41 #include <utility>
42 
43 DetermineTowerBackground::DetermineTowerBackground(const std::string &name)
44  : SubsysReco(name)
45  , _do_flow(0)
46  , _v2(0)
47  , _Psi2(0)
48  , _nStrips(0)
49  , _nTowers(0)
50  , _HCAL_NETA(-1)
51  , _HCAL_NPHI(-1)
52  , _backgroundName("TestTowerBackground")
53  , _seed_type(0)
54  , _seed_jet_D(3.0)
55  , _seed_jet_pt(7.0)
56 {
57  _UE.resize(3, std::vector<float>(1, 0));
58 }
59 
60 int DetermineTowerBackground::InitRun(PHCompositeNode *topNode)
61 {
62  return CreateNode(topNode);
63 }
64 
65 int DetermineTowerBackground::process_event(PHCompositeNode *topNode)
66 {
67  if (Verbosity() > 0)
68  {
69  std::cout << "DetermineTowerBackground::process_event: entering with do_flow = " << _do_flow << ", seed type = " << _seed_type << ", ";
70  if (_seed_type == 0)
71  std::cout << " D = " << _seed_jet_D << std::endl;
72  else if (_seed_type == 1)
73  std::cout << " pT = " << _seed_jet_pt << std::endl;
74  else
75  std::cout << " UNDEFINED seed behavior! " << std::endl;
76  }
77 
78  // clear seed eta/phi positions
79  _seed_eta.resize(0);
80  _seed_phi.resize(0);
81 
82  // pull out the tower containers and geometry objects at the start
83  RawTowerContainer *towersEM3 = nullptr;
84  RawTowerContainer *towersIH3 = nullptr;
85  RawTowerContainer *towersOH3 = nullptr;
86  TowerInfoContainer *towerinfosEM3 = nullptr;
87  TowerInfoContainer *towerinfosIH3 = nullptr;
88  TowerInfoContainer *towerinfosOH3 = nullptr;
89  if (m_use_towerinfo)
90  {
91  towerinfosEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
92  towerinfosIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
93  towerinfosOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
94  }
95  else
96  {
97  towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
98  towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
99  towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
100 
101  if (Verbosity() > 0)
102  {
103  std::cout << "DetermineTowerBackground::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
104  std::cout << "DetermineTowerBackground::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << std::endl;
105  std::cout << "DetermineTowerBackground::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
106  }
107  }
108 
109 
110  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
111  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
112 
113  // seed type 0 is D > 3 R=0.2 jets run on retowerized CEMC
114  if (_seed_type == 0)
115  {
116  JetContainer *reco2_jets;
117  if (m_use_towerinfo)
118  {
119  reco2_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
120  }
121  else
122  {
123  reco2_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_HIRecoSeedsRaw_r02");
124  }
125  if (Verbosity() > 1)
126  std::cout << "DetermineTowerBackground::process_event: examining possible seeds (1st iteration) ... " << std::endl;
127 
128  _index_SeedD = reco2_jets->property_index(Jet::PROPERTY::prop_SeedD);
129  _index_SeedItr = reco2_jets->property_index(Jet::PROPERTY::prop_SeedItr);
130  for (auto this_jet : *reco2_jets) {
131 
132  float this_pt = this_jet->get_pt();
133  float this_phi = this_jet->get_phi();
134  float this_eta = this_jet->get_eta();
135 
136  if (this_jet->get_pt() < 5)
137  {
138  // mark that this jet was not selected as a seed (and did not have D determined)
139  this_jet->set_property(_index_SeedD, 0);
140  this_jet->set_property(_index_SeedItr, 0);
141 
142  continue;
143  }
144 
145  if (Verbosity() > 2)
146  std::cout << "DetermineTowerBackground::process_event: possible seed jet with pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << ", examining constituents..." << std::endl;
147 
148  std::map<int, double> constituent_ETsum;
149 
150  for (const auto& comp : this_jet->get_comp_vec())
151  {
152  int comp_ieta = -1;
153  int comp_iphi = -1;
154  float comp_ET = 0;
155  int comp_T = -99;
156 
157  RawTower *tower;
158  TowerInfo* towerinfo;
159  RawTowerGeom *tower_geom;
160 
161  if (m_use_towerinfo)
162  {
163  if (comp.first == 5 || comp.first == 26)
164  {
165  towerinfo = towerinfosIH3->get_tower_at_channel(comp.second);
166  unsigned int towerkey = towerinfosIH3->encode_key(comp.second);
167  comp_ieta = towerinfosIH3->getTowerEtaBin(towerkey);
168  comp_iphi = towerinfosIH3->getTowerPhiBin(towerkey);
169  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, comp_ieta, comp_iphi);
170  tower_geom = geomIH->get_tower_geometry(key);
171  comp_ET = towerinfo->get_energy() / cosh(tower_geom->get_eta());
172  comp_T = towerinfo->get_time();
173  }
174  else if (comp.first == 7 || comp.first == 27)
175  {
176  towerinfo = towerinfosOH3->get_tower_at_channel(comp.second);
177  unsigned int towerkey = towerinfosOH3->encode_key(comp.second);
178  comp_ieta = towerinfosOH3->getTowerEtaBin(towerkey);
179  comp_iphi = towerinfosOH3->getTowerPhiBin(towerkey);
180  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, comp_ieta, comp_iphi);
181  tower_geom = geomOH->get_tower_geometry(key);
182  comp_ET = towerinfo->get_energy() / cosh(tower_geom->get_eta());
183  comp_T = towerinfo->get_time();
184  }
185  else if (comp.first == 13 || comp.first == 28)
186  {
187  towerinfo = towerinfosEM3->get_tower_at_channel(comp.second);
188  unsigned int towerkey = towerinfosEM3->encode_key(comp.second);
189  comp_ieta = towerinfosEM3->getTowerEtaBin(towerkey);
190  comp_iphi = towerinfosEM3->getTowerPhiBin(towerkey);
191  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, comp_ieta, comp_iphi);
192  tower_geom = geomIH->get_tower_geometry(key);
193  comp_ET = towerinfo->get_energy() / cosh(tower_geom->get_eta());
194  comp_T = towerinfo->get_time();
195  }
196  }
197  else
198  {
199 
200  if (comp.first == 5)
201  {
202  tower = towersIH3->getTower(comp.second);
203  tower_geom = geomIH->get_tower_geometry(tower->get_key());
204 
205  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
206  comp_iphi = geomIH->get_phibin(tower_geom->get_phi());
207  comp_ET = tower->get_energy() / cosh(tower_geom->get_eta());
208  }
209  else if (comp.first == 7)
210  {
211  tower = towersOH3->getTower(comp.second);
212  tower_geom = geomOH->get_tower_geometry(tower->get_key());
213 
214  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
215  comp_iphi = geomIH->get_phibin(tower_geom->get_phi());
216  comp_ET = tower->get_energy() / cosh(tower_geom->get_eta());
217  }
218  else if (comp.first == 13)
219  {
220  tower = towersEM3->getTower(comp.second);
221  tower_geom = geomIH->get_tower_geometry(tower->get_key());
222 
223  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
224  comp_iphi = geomIH->get_phibin(tower_geom->get_phi());
225  comp_ET = tower->get_energy() / cosh(tower_geom->get_eta());
226  }
227  }
228  if(comp_T == -10 || comp_T == -11)
229  {
230  if (Verbosity() > 4)
231  std::cout << "DetermineTowerBackground::process_event: --> --> Skipping constituent in layer " << comp.first << " at ieta / iphi = " << comp_ieta << " / " << comp_iphi << "due to masking" << std::endl;
232  continue;
233  }
234  int comp_ikey = 1000 * comp_ieta + comp_iphi;
235 
236  if (Verbosity() > 4)
237  std::cout << "DetermineTowerBackground::process_event: --> --> constituent in layer " << comp.first << " at ieta / iphi = " << comp_ieta << " / " << comp_iphi << ", filling map with key = " << comp_ikey << " and ET = " << comp_ET << std::endl;
238 
239  constituent_ETsum[comp_ikey] += comp_ET;
240 
241  if (Verbosity() > 4)
242  std::cout << "DetermineTowerBackground::process_event: --> --> ET sum map at key = " << comp_ikey << " now has ET = " << constituent_ETsum[comp_ikey] << std::endl;
243  }
244 
245  // now iterate over constituent_ET sums to find maximum and mean
246  float constituent_max_ET = 0;
247  float constituent_sum_ET = 0;
248  int nconstituents = 0;
249 
250  if (Verbosity() > 4)
251  std::cout << "DetermineTowerBackground::process_event: --> now iterating over map..." << std::endl;
252  for (std::map<int, double>::iterator map_iter = constituent_ETsum.begin(); map_iter != constituent_ETsum.end(); ++map_iter)
253  {
254  if (Verbosity() > 4)
255  std::cout << "DetermineTowerBackground::process_event: --> --> map has key # " << map_iter->first << " and ET = " << map_iter->second << std::endl;
256  nconstituents++;
257  constituent_sum_ET += map_iter->second;
258  if (map_iter->second > constituent_max_ET) constituent_max_ET = map_iter->second;
259  }
260 
261  float mean_constituent_ET = constituent_sum_ET / nconstituents;
262  float seed_D = constituent_max_ET / mean_constituent_ET;
263 
264  // store D value as property for offline analysis / debugging
265  this_jet->set_property(_index_SeedD, seed_D);
266 
267  if (Verbosity() > 3)
268  std::cout << "DetermineTowerBackground::process_event: --> jet has < ET > = " << constituent_sum_ET << " / " << nconstituents << " = " << mean_constituent_ET << ", max-ET = " << constituent_max_ET << ", and D = " << seed_D << std::endl;
269 
270  if (seed_D > _seed_jet_D)
271  {
272  _seed_eta.push_back(this_eta);
273  _seed_phi.push_back(this_phi);
274 
275  // set first iteration seed property
276  this_jet->set_property(_index_SeedItr, 1.0);
277 
278  if (Verbosity() > 1)
279  std::cout << "DetermineTowerBackground::process_event: --> adding seed at eta / phi = " << this_eta << " / " << this_phi << " ( R=0.2 jet with pt = " << this_pt << ", D = " << seed_D << " ) " << std::endl;
280  }
281  else
282  {
283  // mark that this jet was considered but not used as a seed
284  this_jet->set_property(_index_SeedItr, 0.0);
285 
286  if (Verbosity() > 3)
287  std::cout << "DetermineTowerBackground::process_event: --> discarding potential seed at eta / phi = " << this_eta << " / " << this_phi << " ( R=0.2 jet with pt = " << this_pt << ", D = " << seed_D << " ) " << std::endl;
288  }
289  }
290  }
291 
292 
293  // seed type 1 is the set of those jets above which, when their
294  // kinematics are updated for the first background subtraction, have
295  // pT > 20 GeV
296  if (_seed_type == 1)
297  {
298  JetContainer *reco2_jets;
299  if (m_use_towerinfo)
300  {
301  reco2_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
302  }
303  else
304  {
305  reco2_jets = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_HIRecoSeedsSub_r02");
306  }
307  if (Verbosity() > 1)
308  std::cout << "DetermineTowerBackground::process_event: examining possible seeds (2nd iteration) ... " << std::endl;
309 
310  _index_SeedD = reco2_jets->property_index(Jet::PROPERTY::prop_SeedD);
311  _index_SeedItr = reco2_jets->property_index(Jet::PROPERTY::prop_SeedItr);
312  for (auto this_jet : *reco2_jets)
313  {
314  float this_pt = this_jet->get_pt();
315  float this_phi = this_jet->get_phi();
316  float this_eta = this_jet->get_eta();
317 
318  if (this_jet->get_pt() < _seed_jet_pt)
319  {
320  // mark that this jet was considered but not used as a seed
321  this_jet->set_property(_index_SeedItr, 0.0);
322 
323  continue;
324  }
325 
326  _seed_eta.push_back(this_eta);
327  _seed_phi.push_back(this_phi);
328 
329  // set second iteration seed property
330  this_jet->set_property(_index_SeedItr, 2.0);
331 
332  if (Verbosity() > 1)
333  std::cout << "DetermineTowerBackground::process_event: --> adding seed at eta / phi = " << this_eta << " / " << this_phi << " ( R=0.2 jet with pt = " << this_pt << " ) " << std::endl;
334  }
335  }
336 
337  // get the binning from the geometry (different for 1D vs 2D...)
338  if (_HCAL_NETA < 0)
339  {
340  _HCAL_NETA = geomIH->get_etabins();
341  _HCAL_NPHI = geomIH->get_phibins();
342 
343  // resize UE density and energy vectors
344  _UE[0].resize(_HCAL_NETA, 0);
345  _UE[1].resize(_HCAL_NETA, 0);
346  _UE[2].resize(_HCAL_NETA, 0);
347 
348  _EMCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
349  _IHCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
350  _OHCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
351 
352  _EMCAL_T.resize(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
353  _IHCAL_T.resize(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
354  _OHCAL_T.resize(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0));
355 
356  // for flow determination, build up a 1-D phi distribution of
357  // energies from all layers summed together, populated only from eta
358  // strips which do not have any excluded phi towers
359  _FULLCALOFLOW_PHI_E.resize(_HCAL_NPHI, 0);
360  _FULLCALOFLOW_PHI_VAL.resize(_HCAL_NPHI, 0);
361 
362  if (Verbosity() > 0)
363  {
364  std::cout << "DetermineTowerBackground::process_event: setting number of towers in eta / phi: " << _HCAL_NETA << " / " << _HCAL_NPHI << std::endl;
365  }
366  }
367 
368  // reset all maps map
369  for (int ieta = 0; ieta < _HCAL_NETA; ieta++)
370  {
371  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
372  {
373  _EMCAL_E[ieta][iphi] = 0;
374  _IHCAL_E[ieta][iphi] = 0;
375  _OHCAL_E[ieta][iphi] = 0;
376  _EMCAL_T[ieta][iphi] = 0;
377  _IHCAL_T[ieta][iphi] = 0;
378  _OHCAL_T[ieta][iphi] = 0;
379  }
380  }
381 
382  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
383  {
384  _FULLCALOFLOW_PHI_E[iphi] = 0;
385  _FULLCALOFLOW_PHI_VAL[iphi] = 0;
386  }
387 
388 
389 
390  if (m_use_towerinfo)
391  {
392  // iterate over EMCal towerinfos
393  if (!towerinfosEM3 || !towerinfosIH3 || !towerinfosOH3)
394  {
395  std::cout << PHWHERE << "missing tower info object, doing nothing" << std::endl;
396  return Fun4AllReturnCodes::ABORTRUN;
397  }
398  unsigned int nchannels_em = towerinfosEM3->size();
399  for (unsigned int channel = 0; channel < nchannels_em;channel++)
400  {
401  unsigned int key = towerinfosEM3->encode_key(channel);
402  int this_etabin = towerinfosEM3->getTowerEtaBin(key);
403  int this_phibin = towerinfosEM3->getTowerPhiBin(key);
404  float this_E = towerinfosEM3->get_tower_at_channel(channel)->get_energy();
405  int this_T = towerinfosEM3->get_tower_at_channel(channel)->get_time();
406  _EMCAL_E[this_etabin][this_phibin] += this_E;
407  _EMCAL_T[this_etabin][this_phibin] = this_T;
408  }
409 
410  // iterate over IHCal towerinfos
411  unsigned int nchannels_ih = towerinfosIH3->size();
412  for (unsigned int channel = 0; channel < nchannels_ih;channel++)
413  {
414  unsigned int key = towerinfosIH3->encode_key(channel);
415  int this_etabin = towerinfosIH3->getTowerEtaBin(key);
416  int this_phibin = towerinfosIH3->getTowerPhiBin(key);
417  float this_E = towerinfosIH3->get_tower_at_channel(channel)->get_energy();
418  int this_T = towerinfosIH3->get_tower_at_channel(channel)->get_time();
419  _IHCAL_E[this_etabin][this_phibin] += this_E;
420  _IHCAL_T[this_etabin][this_phibin] += this_T;
421  }
422 
423  // iterate over OHCal towerinfos
424  unsigned int nchannels_oh = towerinfosOH3->size();
425  for (unsigned int channel = 0; channel < nchannels_oh;channel++)
426  {
427  unsigned int key = towerinfosOH3->encode_key(channel);
428  int this_etabin = towerinfosOH3->getTowerEtaBin(key);
429  int this_phibin = towerinfosOH3->getTowerPhiBin(key);
430  float this_E = towerinfosOH3->get_tower_at_channel(channel)->get_energy();
431  int this_T = towerinfosOH3->get_tower_at_channel(channel)->get_time();
432  _OHCAL_E[this_etabin][this_phibin] += this_E;
433  _OHCAL_T[this_etabin][this_phibin] += this_T;
434  }
435  }
436  else
437  {
438  // iterate over EMCal towers
439  RawTowerContainer::ConstRange begin_end = towersEM3->getTowers();
440  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
441  {
442  RawTower *tower = rtiter->second;
443 
444  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
445 
446  float this_eta = tower_geom->get_eta();
447  float this_phi = tower_geom->get_phi();
448  int this_etabin = geomIH->get_etabin(this_eta);
449  int this_phibin = geomIH->get_phibin(this_phi);
450  float this_E = tower->get_energy();
451 
452  _EMCAL_E[this_etabin][this_phibin] += this_E;
453 
454  if (Verbosity() > 2 && tower->get_energy() > 1)
455  {
456  std::cout << "DetermineTowerBackground::process_event: EMCal tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << this_etabin << " ) / " << this_phi << " ( " << this_phibin << " ) / " << this_E << std::endl;
457  }
458  }
459 
460  // iterate over IHCal towers
461  RawTowerContainer::ConstRange begin_end_IH = towersIH3->getTowers();
462  for (RawTowerContainer::ConstIterator rtiter = begin_end_IH.first; rtiter != begin_end_IH.second; ++rtiter)
463  {
464  RawTower *tower = rtiter->second;
465  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
466 
467  float this_eta = tower_geom->get_eta();
468  float this_phi = tower_geom->get_phi();
469  int this_etabin = geomIH->get_etabin(this_eta);
470  int this_phibin = geomIH->get_phibin(this_phi);
471  float this_E = tower->get_energy();
472 
473  _IHCAL_E[this_etabin][this_phibin] += this_E;
474 
475  if (Verbosity() > 2 && tower->get_energy() > 1)
476  {
477  std::cout << "DetermineTowerBackground::process_event: IHCal tower at eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << this_etabin << " ) / " << this_phi << " ( " << this_phibin << " ) / " << this_E << std::endl;
478  }
479  }
480 
481  // iterate over OHCal towers
482  RawTowerContainer::ConstRange begin_end_OH = towersOH3->getTowers();
483  for (RawTowerContainer::ConstIterator rtiter = begin_end_OH.first; rtiter != begin_end_OH.second; ++rtiter)
484  {
485  RawTower *tower = rtiter->second;
486  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
487 
488  float this_eta = tower_geom->get_eta();
489  float this_phi = tower_geom->get_phi();
490  int this_etabin = geomOH->get_etabin(this_eta);
491  int this_phibin = geomOH->get_phibin(this_phi);
492  float this_E = tower->get_energy();
493 
494  _OHCAL_E[this_etabin][this_phibin] += this_E;
495 
496  if (Verbosity() > 2 && tower->get_energy() > 1)
497  {
498  std::cout << "DetermineTowerBackground::process_event: OHCal tower at eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << this_etabin << " ) / " << this_phi << " ( " << this_phibin << " ) / " << this_E << std::endl;
499  }
500  }
501  }
502  // first, calculate flow: Psi2 & v2, if enabled
503 
504  _Psi2 = 0;
505  _v2 = 0;
506  _nStrips = 0;
507 
508  if (_do_flow == 0)
509  {
510  if (Verbosity() > 0)
511  {
512  std::cout << "DetermineTowerBackground::process_event: flow not enabled, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) , v2 = " << _v2 << std::endl;
513  }
514  }
515 
516  if (_do_flow >= 1)
517  {
518  // check for the case when every tower is excluded
519  int nStripsAvailableForFlow = 0;
520  int nStripsUnavailableForFlow = 0;
521 
522  for (int layer = 0; layer < 3; layer++)
523  {
524  int local_max_eta = _HCAL_NETA;
525  int local_max_phi = _HCAL_NPHI;
526 
527  for (int eta = 0; eta < local_max_eta; eta++)
528  {
529  bool isAnyTowerExcluded = false;
530 
531  for (int phi = 0; phi < local_max_phi; phi++)
532  {
533  float this_eta = geomIH->get_etacenter(eta);
534  float this_phi = geomIH->get_phicenter(phi);
535 
536  bool isExcluded = false;
537 
538  //if the tower is masked, exclude it
539  int my_T = (layer == 0 ? _EMCAL_T[eta][phi] : (layer == 1 ? _IHCAL_T[eta][phi] : _OHCAL_T[eta][phi]));
540  if(my_T == -10 || my_T == -11) isExcluded = true;
541  for (unsigned int iseed = 0; iseed < _seed_eta.size(); iseed++)
542  {
543  float deta = this_eta - _seed_eta[iseed];
544  float dphi = this_phi - _seed_phi[iseed];
545  if (dphi > M_PI) dphi -= 2 * M_PI;
546  if (dphi < -M_PI) dphi += 2 * M_PI;
547  float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
548  if (dR < 0.4)
549  {
550  isExcluded = true;
551  if (Verbosity() > 10)
552  {
553  float my_E = (layer == 0 ? _EMCAL_E[eta][phi] : (layer == 1 ? _IHCAL_E[eta][phi] : _OHCAL_E[eta][phi]));
554  std::cout << " setting excluded mark for tower with E / eta / phi = " << my_E << " / " << this_eta << " / " << this_phi << " from seed at eta / phi = " << _seed_eta[iseed] << " / " << _seed_phi[iseed] << std::endl;
555  }
556  }
557  }
558 
559  // if even a single tower in this eta strip is excluded, we
560  // can't use the strip for flow determination
561  if (isExcluded)
562  isAnyTowerExcluded = true;
563  } // close phi loop
564 
565  // if this eta strip can be used for flow determination, fill it now
566  if (!isAnyTowerExcluded)
567  {
568  if (Verbosity() > 4)
569  std::cout << " strip at layer " << layer << ", eta " << eta << " has no excluded towers and can be used for flow determination " << std::endl;
570  nStripsAvailableForFlow++;
571 
572  for (int phi = 0; phi < local_max_phi; phi++)
573  {
574  float this_phi = geomIH->get_phicenter(phi);
575 
576  if (layer == 0) _FULLCALOFLOW_PHI_E[phi] += _EMCAL_E[eta][phi];
577  if (layer == 1) _FULLCALOFLOW_PHI_E[phi] += _IHCAL_E[eta][phi];
578  if (layer == 2) _FULLCALOFLOW_PHI_E[phi] += _OHCAL_E[eta][phi];
579 
580  _FULLCALOFLOW_PHI_VAL[phi] = this_phi; // should really set this globally only one time
581  }
582  }
583  else
584  {
585  if (Verbosity() > 4)
586  std::cout << " strip at layer " << layer << ", eta " << eta << " DOES have excluded towers and CANNOT be used for flow determination " << std::endl;
587  nStripsUnavailableForFlow++;
588  }
589 
590  } // close eta loop
591  } // close layer loop
592 
593  // flow determination
594 
595  float Q_x = 0;
596  float Q_y = 0;
597  float E = 0;
598 
599  if (Verbosity() > 0)
600  std::cout << "DetermineTowerBackground::process_event: # of strips (summed over layers) available / unavailable for flow determination: " << nStripsAvailableForFlow << " / " << nStripsUnavailableForFlow << std::endl;
601 
602  if (nStripsAvailableForFlow > 0)
603  {
604  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
605  {
606  E += _FULLCALOFLOW_PHI_E[iphi];
607  Q_x += _FULLCALOFLOW_PHI_E[iphi] * cos(2 * _FULLCALOFLOW_PHI_VAL[iphi]);
608  Q_y += _FULLCALOFLOW_PHI_E[iphi] * sin(2 * _FULLCALOFLOW_PHI_VAL[iphi]);
609  }
610 
611  if (_do_flow == 1)
612  {
613  _Psi2 = atan2(Q_y, Q_x) / 2.0;
614  }
615  else if (_do_flow == 2)
616  {
617  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
618 
619  if (!truthinfo)
620  {
621  std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl;
622  return -1;
623  }
624 
625  PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
626 
627  float Hijing_Qx = 0, Hijing_Qy = 0;
628 
629  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
630  {
631  PHG4Particle *g4particle = iter->second;
632 
633  if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0) continue;
634 
635  TLorentzVector t;
636  t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e());
637 
638  float truth_pt = t.Pt();
639  if (truth_pt < 0.4) continue;
640  float truth_eta = t.Eta();
641  if (fabs(truth_eta) > 1.1) continue;
642  float truth_phi = t.Phi();
643  int truth_pid = g4particle->get_pid();
644 
645  if (Verbosity() > 10)
646  std::cout << "DetermineTowerBackground::process_event: determining truth flow, using particle w/ pt / eta / phi " << truth_pt << " / " << truth_eta << " / " << truth_phi << " , embed / PID = " << truthinfo->isEmbeded(g4particle->get_track_id()) << " / " << truth_pid << std::endl;
647 
648  Hijing_Qx += truth_pt * cos(2 * truth_phi);
649  Hijing_Qy += truth_pt * sin(2 * truth_phi);
650  }
651 
652  _Psi2 = atan2(Hijing_Qy, Hijing_Qx) / 2.0;
653 
654  if (Verbosity() > 0)
655  std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl;
656  }
657 
658  // determine v2 from calo regardless of origin of Psi2
659  double sum_cos2dphi = 0;
660  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
661  {
662  sum_cos2dphi += _FULLCALOFLOW_PHI_E[iphi] * cos(2 * (_FULLCALOFLOW_PHI_VAL[iphi] - _Psi2));
663  }
664 
665  _v2 = sum_cos2dphi / E;
666 
667  _nStrips = nStripsAvailableForFlow;
668  }
669  else
670  {
671  _Psi2 = 0;
672  _v2 = 0;
673  _nStrips = 0;
674  if (Verbosity() > 0)
675  std::cout << "DetermineTowerBackground::process_event: no full strips available for flow modulation, setting v2 and Psi = 0" << std::endl;
676  }
677 
678  if (Verbosity() > 0)
679  {
680  std::cout << "DetermineTowerBackground::process_event: unnormalized Q vector (Qx, Qy) = ( " << Q_x << ", " << Q_y << " ) with Sum E_i = " << E << std::endl;
681  std::cout << "DetermineTowerBackground::process_event: Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi " << (_do_flow == 2 ? "from Hijing " : "") << ") , v2 = " << _v2 << " ( using " << _nStrips << " ) " << std::endl;
682  }
683  } // if do flow
684 
685  // now calculate energy densities...
686  _nTowers = 0; // store how many towers were used to determine bkg
687 
688  // starting with the EMCal first...
689  for (int layer = 0; layer < 3; layer++)
690  {
691  int local_max_eta = _HCAL_NETA;
692  int local_max_phi = _HCAL_NPHI;
693 
694  for (int eta = 0; eta < local_max_eta; eta++)
695  {
696  float total_E = 0;
697  int total_tower = 0;
698 
699  for (int phi = 0; phi < local_max_phi; phi++)
700  {
701  float this_eta = geomIH->get_etacenter(eta);
702  float this_phi = geomIH->get_phicenter(phi);
703 
704  bool isExcluded = false;
705 
706  float my_E = (layer == 0 ? _EMCAL_E[eta][phi] : (layer == 1 ? _IHCAL_E[eta][phi] : _OHCAL_E[eta][phi]));
707  int my_T = (layer == 0 ? _EMCAL_T[eta][phi] : (layer == 1 ? _IHCAL_T[eta][phi] : _OHCAL_T[eta][phi]));
708  //if the tower is masked (energy identically zero), exclude it
709  if(my_T == -10 || my_T == -11)
710  {
711  isExcluded = true;
712  if (Verbosity() > 10) std::cout << " tower in layer " << layer << " at eta / phi = " << this_eta << " / " << this_phi << " with E = " << my_E << " excluded due to masking" << std::endl;
713  }
714  for (unsigned int iseed = 0; iseed < _seed_eta.size(); iseed++)
715  {
716  float deta = this_eta - _seed_eta[iseed];
717  float dphi = this_phi - _seed_phi[iseed];
718  if (dphi > M_PI) dphi -= 2 * M_PI;
719  if (dphi < -M_PI) dphi += 2 * M_PI;
720  float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
721  if (dR < 0.4)
722  {
723  isExcluded = true;
724  if (Verbosity() > 10) std::cout << " setting excluded mark from seed at eta / phi = " << _seed_eta[iseed] << " / " << _seed_phi[iseed] << std::endl;
725  }
726  }
727 
728  if (!isExcluded)
729  {
730  if (layer == 0) total_E += _EMCAL_E[eta][phi] / (1 + 2 * _v2 * cos(2 * (this_phi - _Psi2)));
731  if (layer == 1) total_E += _IHCAL_E[eta][phi] / (1 + 2 * _v2 * cos(2 * (this_phi - _Psi2)));
732  if (layer == 2) total_E += _OHCAL_E[eta][phi] / (1 + 2 * _v2 * cos(2 * (this_phi - _Psi2)));
733  total_tower++; // towers in this eta range & layer
734  _nTowers++; // towers in entire calorimeter
735  }
736  else
737  {
738  if (Verbosity() > 10) std::cout << " tower at eta / phi = " << this_eta << " / " << this_phi << " with E = " << total_E << " excluded due to seed " << std::endl;
739  }
740  }
741 
742  std::pair<float, float> etabounds = geomIH->get_etabounds(eta);
743  std::pair<float, float> phibounds = geomIH->get_phibounds(0);
744 
745  float deta = etabounds.second - etabounds.first;
746  float dphi = phibounds.second - phibounds.first;
747  float total_area = total_tower * deta * dphi;
748 
749  _UE[layer].at(eta) = total_E / total_tower;
750 
751  if (Verbosity() > 3)
752  {
753  std::cout << "DetermineTowerBackground::process_event: at layer / eta index ( eta range ) = " << layer << " / " << eta << " ( " << etabounds.first << " - " << etabounds.second << " ) , total E / total Ntower / total area = " << total_E << " / " << total_tower << " / " << total_area << " , UE per tower = " << total_E / total_tower << std::endl;
754  }
755  }
756  }
757 
758  if (Verbosity() > 0)
759  {
760  for (int layer = 0; layer < 3; layer++)
761  {
762  std::cout << "DetermineTowerBackground::process_event: summary UE in layer " << layer << " : ";
763  for (int eta = 0; eta < _HCAL_NETA; eta++) std::cout << _UE[layer].at(eta) << " , ";
764  std::cout << std::endl;
765  }
766  }
767 
768  //
769 
770  FillNode(topNode);
771 
772  if (Verbosity() > 0) std::cout << "DetermineTowerBackground::process_event: exiting" << std::endl;
773 
774  return Fun4AllReturnCodes::EVENT_OK;
775 }
776 
777 int DetermineTowerBackground::CreateNode(PHCompositeNode *topNode)
778 {
779  PHNodeIterator iter(topNode);
780 
781  // Looking for the DST node
782  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
783  if (!dstNode)
784  {
785  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
786  return Fun4AllReturnCodes::ABORTRUN;
787  }
788 
789  // store the jet background stuff under a sub-node directory
790  PHCompositeNode *bkgNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "JETBACKGROUND"));
791  if (!bkgNode)
792  {
793  bkgNode = new PHCompositeNode("JETBACKGROUND");
794  dstNode->addNode(bkgNode);
795  }
796 
797  // create the TowerBackground node...
798  TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, _backgroundName);
799  if (!towerbackground)
800  {
801  towerbackground = new TowerBackgroundv1();
802  PHIODataNode<PHObject> *bkgDataNode = new PHIODataNode<PHObject>(towerbackground, _backgroundName, "PHObject");
803  bkgNode->addNode(bkgDataNode);
804  }
805  else
806  {
807  std::cout << PHWHERE << "::ERROR - " << _backgroundName << " pre-exists, but should not" << std::endl;
808  exit(-1);
809  }
810 
811  return Fun4AllReturnCodes::EVENT_OK;
812 }
813 
814 void DetermineTowerBackground::FillNode(PHCompositeNode *topNode)
815 {
816  TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, _backgroundName);
817  if (!towerbackground)
818  {
819  std::cout << " ERROR -- can't find TowerBackground node after it should have been created" << std::endl;
820  return;
821  }
822  else
823  {
824  towerbackground->set_UE(0, _UE[0]);
825  towerbackground->set_UE(1, _UE[1]);
826  towerbackground->set_UE(2, _UE[2]);
827 
828  towerbackground->set_v2(_v2);
829 
830  towerbackground->set_Psi2(_Psi2);
831 
832  towerbackground->set_nStripsUsedForFlow(_nStrips);
833 
834  towerbackground->set_nTowersUsedForBkg(_nTowers);
835  }
836 
837  return;
838 }