ECCE @ EIC 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
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 <g4jets/Jet.h>
12 #include <g4jets/JetMap.h>
13 
14 #include <g4main/PHG4Particle.h>
16 
18 #include <fun4all/SubsysReco.h>
19 
20 #include <phool/PHCompositeNode.h>
21 #include <phool/PHIODataNode.h>
22 #include <phool/PHNode.h>
23 #include <phool/PHNodeIterator.h>
24 #include <phool/PHObject.h>
25 #include <phool/getClass.h>
26 #include <phool/phool.h>
27 
28 #include <TLorentzVector.h>
29 
30 // standard includes
31 #include <cmath>
32 #include <cstdlib>
33 #include <iomanip>
34 #include <iostream>
35 #include <map>
36 #include <memory>
37 #include <vector>
38 #include <utility>
39 
41  : SubsysReco(name)
42  , _do_flow(0)
43  , _v2(0)
44  , _Psi2(0)
45  , _nStrips(0)
46  , _nTowers(0)
47  , _HCAL_NETA(-1)
48  , _HCAL_NPHI(-1)
49  , _backgroundName("TestTowerBackground")
50  , _seed_type(0)
51  , _seed_jet_D(3.0)
52  , _seed_jet_pt(7.0)
53 {
54  _UE.resize(3, std::vector<float>(1, 0));
55 }
56 
58 {
59  return CreateNode(topNode);
60 }
61 
63 {
64  if (Verbosity() > 0)
65  {
66  std::cout << "DetermineTowerBackground::process_event: entering with do_flow = " << _do_flow << ", seed type = " << _seed_type << ", ";
67  if (_seed_type == 0)
68  std::cout << " D = " << _seed_jet_D << std::endl;
69  else if (_seed_type == 1)
70  std::cout << " pT = " << _seed_jet_pt << std::endl;
71  else
72  std::cout << " UNDEFINED seed behavior! " << std::endl;
73  }
74 
75  // clear seed eta/phi positions
76  _seed_eta.resize(0);
77  _seed_phi.resize(0);
78 
79  // pull out the tower containers and geometry objects at the start
80  RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
81  RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
82  RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
83  if (Verbosity() > 0)
84  {
85  std::cout << "DetermineTowerBackground::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
86  std::cout << "DetermineTowerBackground::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << std::endl;
87  std::cout << "DetermineTowerBackground::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
88  }
89 
90  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
91  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
92 
93  // seed type 0 is D > 3 R=0.2 jets run on retowerized CEMC
94  if (_seed_type == 0)
95  {
96  JetMap *reco2_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_HIRecoSeedsRaw_r02");
97 
98  if (Verbosity() > 1)
99  std::cout << "DetermineTowerBackground::process_event: examining possible seeds (1st iteration) ... " << std::endl;
100 
101  for (JetMap::Iter iter = reco2_jets->begin(); iter != reco2_jets->end(); ++iter)
102  {
103  Jet *this_jet = iter->second;
104 
105  float this_pt = this_jet->get_pt();
106  float this_phi = this_jet->get_phi();
107  float this_eta = this_jet->get_eta();
108 
109  if (this_jet->get_pt() < 5)
110  {
111  // mark that this jet was not selected as a seed (and did not have D determined)
112  this_jet->set_property(Jet::PROPERTY::prop_SeedD, 0);
113  this_jet->set_property(Jet::PROPERTY::prop_SeedItr, 0);
114 
115  continue;
116  }
117 
118  if (Verbosity() > 2)
119  std::cout << "DetermineTowerBackground::process_event: possible seed jet with pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << ", examining constituents..." << std::endl;
120 
121  std::map<int, double> constituent_ETsum;
122 
123  for (Jet::ConstIter comp = this_jet->begin_comp(); comp != this_jet->end_comp(); ++comp)
124  {
125  int comp_ieta = -1;
126  int comp_iphi = -1;
127  float comp_ET = 0;
128 
129  RawTower *tower;
130  RawTowerGeom *tower_geom;
131 
132  if ((*comp).first == 5)
133  {
134  tower = towersIH3->getTower((*comp).second);
135  tower_geom = geomIH->get_tower_geometry(tower->get_key());
136 
137  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
138  comp_iphi = geomIH->get_phibin(tower_geom->get_phi());
139  comp_ET = tower->get_energy() / cosh(tower_geom->get_eta());
140  }
141  else if ((*comp).first == 7)
142  {
143  tower = towersOH3->getTower((*comp).second);
144  tower_geom = geomOH->get_tower_geometry(tower->get_key());
145 
146  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
147  comp_iphi = geomIH->get_phibin(tower_geom->get_phi());
148  comp_ET = tower->get_energy() / cosh(tower_geom->get_eta());
149  }
150  else if ((*comp).first == 13)
151  {
152  tower = towersEM3->getTower((*comp).second);
153  tower_geom = geomIH->get_tower_geometry(tower->get_key());
154 
155  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
156  comp_iphi = geomIH->get_phibin(tower_geom->get_phi());
157  comp_ET = tower->get_energy() / cosh(tower_geom->get_eta());
158  }
159 
160  int comp_ikey = 1000 * comp_ieta + comp_iphi;
161 
162  if (Verbosity() > 4)
163  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;
164 
165  constituent_ETsum[comp_ikey] += comp_ET;
166 
167  if (Verbosity() > 4)
168  std::cout << "DetermineTowerBackground::process_event: --> --> ET sum map at key = " << comp_ikey << " now has ET = " << constituent_ETsum[comp_ikey] << std::endl;
169  }
170 
171  // now iterate over constituent_ET sums to find maximum and mean
172  float constituent_max_ET = 0;
173  float constituent_sum_ET = 0;
174  int nconstituents = 0;
175 
176  if (Verbosity() > 4)
177  std::cout << "DetermineTowerBackground::process_event: --> now iterating over map..." << std::endl;
178  for (std::map<int, double>::iterator map_iter = constituent_ETsum.begin(); map_iter != constituent_ETsum.end(); ++map_iter)
179  {
180  if (Verbosity() > 4)
181  std::cout << "DetermineTowerBackground::process_event: --> --> map has key # " << map_iter->first << " and ET = " << map_iter->second << std::endl;
182  nconstituents++;
183  constituent_sum_ET += map_iter->second;
184  if (map_iter->second > constituent_max_ET) constituent_max_ET = map_iter->second;
185  }
186 
187  float mean_constituent_ET = constituent_sum_ET / nconstituents;
188  float seed_D = constituent_max_ET / mean_constituent_ET;
189 
190  // store D value as property for offline analysis / debugging
191  this_jet->set_property(Jet::PROPERTY::prop_SeedD, seed_D);
192 
193  if (Verbosity() > 3)
194  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;
195 
196  if (seed_D > _seed_jet_D)
197  {
198  _seed_eta.push_back(this_eta);
199  _seed_phi.push_back(this_phi);
200 
201  // set first iteration seed property
202  this_jet->set_property(Jet::PROPERTY::prop_SeedItr, 1.0);
203 
204  if (Verbosity() > 1)
205  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;
206  }
207  else
208  {
209  // mark that this jet was considered but not used as a seed
210  this_jet->set_property(Jet::PROPERTY::prop_SeedItr, 0.0);
211 
212  if (Verbosity() > 3)
213  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;
214  }
215  }
216  }
217 
218  // seed type 1 is the set of those jets above which, when their
219  // kinematics are updated for the first background subtraction, have
220  // pT > 20 GeV
221  if (_seed_type == 1)
222  {
223  JetMap *reco2_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_HIRecoSeedsSub_r02");
224 
225  if (Verbosity() > 1)
226  std::cout << "DetermineTowerBackground::process_event: examining possible seeds (2nd iteration) ... " << std::endl;
227 
228  for (JetMap::Iter iter = reco2_jets->begin(); iter != reco2_jets->end(); ++iter)
229  {
230  Jet *this_jet = iter->second;
231 
232  float this_pt = this_jet->get_pt();
233  float this_phi = this_jet->get_phi();
234  float this_eta = this_jet->get_eta();
235 
236  if (this_jet->get_pt() < _seed_jet_pt)
237  {
238  // mark that this jet was considered but not used as a seed
239  this_jet->set_property(Jet::PROPERTY::prop_SeedItr, 0.0);
240 
241  continue;
242  }
243 
244  _seed_eta.push_back(this_eta);
245  _seed_phi.push_back(this_phi);
246 
247  // set second iteration seed property
248  this_jet->set_property(Jet::PROPERTY::prop_SeedItr, 2.0);
249 
250  if (Verbosity() > 1)
251  std::cout << "DetermineTowerBackground::process_event: --> adding seed at eta / phi = " << this_eta << " / " << this_phi << " ( R=0.2 jet with pt = " << this_pt << " ) " << std::endl;
252  }
253  }
254 
255  // get the binning from the geometry (different for 1D vs 2D...)
256  if (_HCAL_NETA < 0)
257  {
258  _HCAL_NETA = geomIH->get_etabins();
259  _HCAL_NPHI = geomIH->get_phibins();
260 
261  // resize UE density and energy vectors
262  _UE[0].resize(_HCAL_NETA, 0);
263  _UE[1].resize(_HCAL_NETA, 0);
264  _UE[2].resize(_HCAL_NETA, 0);
265 
266  _EMCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
267  _IHCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
268  _OHCAL_E.resize(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0));
269 
270  // for flow determination, build up a 1-D phi distribution of
271  // energies from all layers summed together, populated only from eta
272  // strips which do not have any excluded phi towers
273  _FULLCALOFLOW_PHI_E.resize(_HCAL_NPHI, 0);
275 
276  if (Verbosity() > 0)
277  {
278  std::cout << "DetermineTowerBackground::process_event: setting number of towers in eta / phi: " << _HCAL_NETA << " / " << _HCAL_NPHI << std::endl;
279  }
280  }
281 
282  // reset all maps map
283  for (int ieta = 0; ieta < _HCAL_NETA; ieta++)
284  {
285  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
286  {
287  _EMCAL_E[ieta][iphi] = 0;
288  _IHCAL_E[ieta][iphi] = 0;
289  _OHCAL_E[ieta][iphi] = 0;
290  }
291  }
292 
293  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
294  {
295  _FULLCALOFLOW_PHI_E[iphi] = 0;
296  _FULLCALOFLOW_PHI_VAL[iphi] = 0;
297  }
298 
299  // iterate over EMCal towers
300  RawTowerContainer::ConstRange begin_end = towersEM3->getTowers();
301  for (RawTowerContainer::ConstIterator rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
302  {
303  RawTower *tower = rtiter->second;
304 
305  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
306 
307  float this_eta = tower_geom->get_eta();
308  float this_phi = tower_geom->get_phi();
309  int this_etabin = geomIH->get_etabin(this_eta);
310  int this_phibin = geomIH->get_phibin(this_phi);
311  float this_E = tower->get_energy();
312 
313  _EMCAL_E[this_etabin][this_phibin] += this_E;
314 
315  if (Verbosity() > 2 && tower->get_energy() > 1)
316  {
317  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;
318  }
319  }
320 
321  // iterate over IHCal towers
322  RawTowerContainer::ConstRange begin_end_IH = towersIH3->getTowers();
323  for (RawTowerContainer::ConstIterator rtiter = begin_end_IH.first; rtiter != begin_end_IH.second; ++rtiter)
324  {
325  RawTower *tower = rtiter->second;
326  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
327 
328  float this_eta = tower_geom->get_eta();
329  float this_phi = tower_geom->get_phi();
330  int this_etabin = geomIH->get_etabin(this_eta);
331  int this_phibin = geomIH->get_phibin(this_phi);
332  float this_E = tower->get_energy();
333 
334  _IHCAL_E[this_etabin][this_phibin] += this_E;
335 
336  if (Verbosity() > 2 && tower->get_energy() > 1)
337  {
338  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;
339  }
340  }
341 
342  // iterate over OHCal towers
343  RawTowerContainer::ConstRange begin_end_OH = towersOH3->getTowers();
344  for (RawTowerContainer::ConstIterator rtiter = begin_end_OH.first; rtiter != begin_end_OH.second; ++rtiter)
345  {
346  RawTower *tower = rtiter->second;
347  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
348 
349  float this_eta = tower_geom->get_eta();
350  float this_phi = tower_geom->get_phi();
351  int this_etabin = geomOH->get_etabin(this_eta);
352  int this_phibin = geomOH->get_phibin(this_phi);
353  float this_E = tower->get_energy();
354 
355  _OHCAL_E[this_etabin][this_phibin] += this_E;
356 
357  if (Verbosity() > 2 && tower->get_energy() > 1)
358  {
359  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;
360  }
361  }
362 
363  // first, calculate flow: Psi2 & v2, if enabled
364 
365  _Psi2 = 0;
366  _v2 = 0;
367  _nStrips = 0;
368 
369  if (_do_flow == 0)
370  {
371  if (Verbosity() > 0)
372  {
373  std::cout << "DetermineTowerBackground::process_event: flow not enabled, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / 3.14159 << " * pi ) , v2 = " << _v2 << std::endl;
374  }
375  }
376 
377  if (_do_flow >= 1)
378  {
379  // check for the case when every tower is excluded
380  int nStripsAvailableForFlow = 0;
381  int nStripsUnavailableForFlow = 0;
382 
383  for (int layer = 0; layer < 3; layer++)
384  {
385  int local_max_eta = _HCAL_NETA;
386  int local_max_phi = _HCAL_NPHI;
387 
388  for (int eta = 0; eta < local_max_eta; eta++)
389  {
390  bool isAnyTowerExcluded = false;
391 
392  for (int phi = 0; phi < local_max_phi; phi++)
393  {
394  float this_eta = geomIH->get_etacenter(eta);
395  float this_phi = geomIH->get_phicenter(phi);
396 
397  bool isExcluded = false;
398 
399  for (unsigned int iseed = 0; iseed < _seed_eta.size(); iseed++)
400  {
401  float deta = this_eta - _seed_eta[iseed];
402  float dphi = this_phi - _seed_phi[iseed];
403  if (dphi > 3.14159) dphi -= 2 * 3.14159;
404  if (dphi < -3.14159) dphi += 2 * 3.14159;
405  float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
406  if (dR < 0.4)
407  {
408  isExcluded = true;
409  if (Verbosity() > 10)
410  {
411  float my_E = (layer == 0 ? _EMCAL_E[eta][phi] : (layer == 1 ? _IHCAL_E[eta][phi] : _OHCAL_E[eta][phi]));
412 
413  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;
414  }
415  }
416  }
417 
418  // if even a single tower in this eta strip is excluded, we
419  // can't use the strip for flow determination
420  if (isExcluded)
421  isAnyTowerExcluded = true;
422  } // close phi loop
423 
424  // if this eta strip can be used for flow determination, fill it now
425  if (!isAnyTowerExcluded)
426  {
427  if (Verbosity() > 4)
428  std::cout << " strip at layer " << layer << ", eta " << eta << " has no excluded towers and can be used for flow determination " << std::endl;
429  nStripsAvailableForFlow++;
430 
431  for (int phi = 0; phi < local_max_phi; phi++)
432  {
433  float this_phi = geomIH->get_phicenter(phi);
434 
435  if (layer == 0) _FULLCALOFLOW_PHI_E[phi] += _EMCAL_E[eta][phi];
436  if (layer == 1) _FULLCALOFLOW_PHI_E[phi] += _IHCAL_E[eta][phi];
437  if (layer == 2) _FULLCALOFLOW_PHI_E[phi] += _OHCAL_E[eta][phi];
438 
439  _FULLCALOFLOW_PHI_VAL[phi] = this_phi; // should really set this globally only one time
440  }
441  }
442  else
443  {
444  if (Verbosity() > 4)
445  std::cout << " strip at layer " << layer << ", eta " << eta << " DOES have excluded towers and CANNOT be used for flow determination " << std::endl;
446  nStripsUnavailableForFlow++;
447  }
448 
449  } // close eta loop
450  } // close layer loop
451 
452  // flow determination
453 
454  float Q_x = 0;
455  float Q_y = 0;
456  float E = 0;
457 
458  if (Verbosity() > 0)
459  std::cout << "DetermineTowerBackground::process_event: # of strips (summed over layers) available / unavailable for flow determination: " << nStripsAvailableForFlow << " / " << nStripsUnavailableForFlow << std::endl;
460 
461  if (nStripsAvailableForFlow > 0)
462  {
463  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
464  {
465  E += _FULLCALOFLOW_PHI_E[iphi];
466  Q_x += _FULLCALOFLOW_PHI_E[iphi] * cos(2 * _FULLCALOFLOW_PHI_VAL[iphi]);
467  Q_y += _FULLCALOFLOW_PHI_E[iphi] * sin(2 * _FULLCALOFLOW_PHI_VAL[iphi]);
468  }
469 
470  if (_do_flow == 1)
471  {
472  _Psi2 = atan2(Q_y, Q_x) / 2.0;
473  }
474  else if (_do_flow == 2)
475  {
476  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
477 
478  if (!truthinfo)
479  {
480  std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl;
481  return -1;
482  }
483 
485 
486  float Hijing_Qx = 0, Hijing_Qy = 0;
487 
488  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
489  {
490  PHG4Particle *g4particle = iter->second;
491 
492  if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0) continue;
493 
494  TLorentzVector t;
495  t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e());
496 
497  float truth_pt = t.Pt();
498  if (truth_pt < 0.4) continue;
499  float truth_eta = t.Eta();
500  if (fabs(truth_eta) > 1.1) continue;
501  float truth_phi = t.Phi();
502  int truth_pid = g4particle->get_pid();
503 
504  if (Verbosity() > 10)
505  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;
506 
507  Hijing_Qx += truth_pt * cos(2 * truth_phi);
508  Hijing_Qy += truth_pt * sin(2 * truth_phi);
509  }
510 
511  _Psi2 = atan2(Hijing_Qy, Hijing_Qx) / 2.0;
512 
513  if (Verbosity() > 0)
514  std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / 3.14159 << " * pi ) " << std::endl;
515  }
516 
517  // determine v2 from calo regardless of origin of Psi2
518  double sum_cos2dphi = 0;
519  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
520  {
521  sum_cos2dphi += _FULLCALOFLOW_PHI_E[iphi] * cos(2 * (_FULLCALOFLOW_PHI_VAL[iphi] - _Psi2));
522  }
523 
524  _v2 = sum_cos2dphi / E;
525 
526  _nStrips = nStripsAvailableForFlow;
527  }
528  else
529  {
530  _Psi2 = 0;
531  _v2 = 0;
532  _nStrips = 0;
533  if (Verbosity() > 0)
534  std::cout << "DetermineTowerBackground::process_event: no full strips available for flow modulation, setting v2 and Psi = 0" << std::endl;
535  }
536 
537  if (Verbosity() > 0)
538  {
539  std::cout << "DetermineTowerBackground::process_event: unnormalized Q vector (Qx, Qy) = ( " << Q_x << ", " << Q_y << " ) with Sum E_i = " << E << std::endl;
540  std::cout << "DetermineTowerBackground::process_event: Psi2 = " << _Psi2 << " ( " << _Psi2 / 3.14159 << " * pi " << (_do_flow == 2 ? "from Hijing " : "") << ") , v2 = " << _v2 << " ( using " << _nStrips << " ) " << std::endl;
541  }
542  } // if do flow
543 
544  // now calculate energy densities...
545  _nTowers = 0; // store how many towers were used to determine bkg
546 
547  // starting with the EMCal first...
548  for (int layer = 0; layer < 3; layer++)
549  {
550  int local_max_eta = _HCAL_NETA;
551  int local_max_phi = _HCAL_NPHI;
552 
553  for (int eta = 0; eta < local_max_eta; eta++)
554  {
555  float total_E = 0;
556  int total_tower = 0;
557 
558  for (int phi = 0; phi < local_max_phi; phi++)
559  {
560  float this_eta = geomIH->get_etacenter(eta);
561  float this_phi = geomIH->get_phicenter(phi);
562 
563  bool isExcluded = false;
564 
565  for (unsigned int iseed = 0; iseed < _seed_eta.size(); iseed++)
566  {
567  float deta = this_eta - _seed_eta[iseed];
568  float dphi = this_phi - _seed_phi[iseed];
569  if (dphi > 3.14159) dphi -= 2 * 3.14159;
570  if (dphi < -3.14159) dphi += 2 * 3.14159;
571  float dR = sqrt(pow(deta, 2) + pow(dphi, 2));
572  if (dR < 0.4)
573  {
574  isExcluded = true;
575  if (Verbosity() > 10) std::cout << " setting excluded mark from seed at eta / phi = " << _seed_eta[iseed] << " / " << _seed_phi[iseed] << std::endl;
576  }
577  }
578 
579  if (!isExcluded)
580  {
581  if (layer == 0) total_E += _EMCAL_E[eta][phi] / (1 + 2 * _v2 * cos(2 * (this_phi - _Psi2)));
582  if (layer == 1) total_E += _IHCAL_E[eta][phi] / (1 + 2 * _v2 * cos(2 * (this_phi - _Psi2)));
583  if (layer == 2) total_E += _OHCAL_E[eta][phi] / (1 + 2 * _v2 * cos(2 * (this_phi - _Psi2)));
584  total_tower++; // towers in this eta range & layer
585  _nTowers++; // towers in entire calorimeter
586  }
587  else
588  {
589  if (Verbosity() > 10) std::cout << " tower at eta / phi = " << this_eta << " / " << this_phi << " with E = " << total_E << " excluded due to seed " << std::endl;
590  }
591  }
592 
593  std::pair<float, float> etabounds = geomIH->get_etabounds(eta);
594  std::pair<float, float> phibounds = geomIH->get_phibounds(0);
595 
596  float deta = etabounds.second - etabounds.first;
597  float dphi = phibounds.second - phibounds.first;
598  float total_area = total_tower * deta * dphi;
599  _UE[layer].at(eta) = total_E / total_tower;
600 
601  if (Verbosity() > 3)
602  {
603  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;
604  }
605  }
606  }
607 
608  if (Verbosity() > 0)
609  {
610  for (int layer = 0; layer < 3; layer++)
611  {
612  std::cout << "DetermineTowerBackground::process_event: summary UE in layer " << layer << " : ";
613  for (int eta = 0; eta < _HCAL_NETA; eta++) std::cout << _UE[layer].at(eta) << " , ";
614  std::cout << std::endl;
615  }
616  }
617 
618  //
619 
620  FillNode(topNode);
621 
622  if (Verbosity() > 0) std::cout << "DetermineTowerBackground::process_event: exiting" << std::endl;
623 
625 }
626 
628 {
629  PHNodeIterator iter(topNode);
630 
631  // Looking for the DST node
632  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
633  if (!dstNode)
634  {
635  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
637  }
638 
639  // store the jet background stuff under a sub-node directory
640  PHCompositeNode *bkgNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "JETBACKGROUND"));
641  if (!bkgNode)
642  {
643  bkgNode = new PHCompositeNode("JETBACKGROUND");
644  dstNode->addNode(bkgNode);
645  }
646 
647  // create the TowerBackground node...
648  TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, _backgroundName);
649  if (!towerbackground)
650  {
651  towerbackground = new TowerBackgroundv1();
652  PHIODataNode<PHObject> *bkgDataNode = new PHIODataNode<PHObject>(towerbackground, _backgroundName, "PHObject");
653  bkgNode->addNode(bkgDataNode);
654  }
655  else
656  {
657  std::cout << PHWHERE << "::ERROR - " << _backgroundName << " pre-exists, but should not" << std::endl;
658  exit(-1);
659  }
660 
662 }
663 
665 {
666  TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, _backgroundName);
667  if (!towerbackground)
668  {
669  std::cout << " ERROR -- can't find TowerBackground node after it should have been created" << std::endl;
670  return;
671  }
672  else
673  {
674  towerbackground->set_UE(0, _UE[0]);
675  towerbackground->set_UE(1, _UE[1]);
676  towerbackground->set_UE(2, _UE[2]);
677 
678  towerbackground->set_v2(_v2);
679 
680  towerbackground->set_Psi2(_Psi2);
681 
682  towerbackground->set_nStripsUsedForFlow(_nStrips);
683 
684  towerbackground->set_nTowersUsedForBkg(_nTowers);
685  }
686 
687  return;
688 }