EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SubtractTowersCS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SubtractTowersCS.cc
1 #include "SubtractTowersCS.h"
2 
3 #include "TowerBackground.h"
4 
5 #include <calobase/RawTower.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerDefs.h>
8 #include <calobase/RawTowerGeom.h>
9 #include <calobase/RawTowerGeomContainer.h>
10 #include <calobase/RawTowerv1.h>
11 
12 #include <fun4all/Fun4AllReturnCodes.h>
13 #include <fun4all/SubsysReco.h>
14 
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHIODataNode.h>
17 #include <phool/PHNode.h>
18 #include <phool/PHNodeIterator.h>
19 #include <phool/PHObject.h>
20 #include <phool/getClass.h>
21 #include <phool/phool.h>
22 
23 #include <fastjet/PseudoJet.hh>
24 #include <fastjet/contrib/ConstituentSubtractor.hh>
25 
26 // standard includes
27 #include <cmath>
28 #include <iostream>
29 #include <map>
30 #include <memory>
31 #include <utility>
32 #include <vector>
33 
34 SubtractTowersCS::SubtractTowersCS(const std::string &name)
35  : SubsysReco(name)
36  , _use_flow_modulation(false)
37  , _alpha(1)
38  , _DeltaRmax(0.3)
39 {
40 }
41 
42 int SubtractTowersCS::InitRun(PHCompositeNode *topNode)
43 {
44  CreateNode(topNode);
45 
46  return Fun4AllReturnCodes::EVENT_OK;
47 }
48 
49 int SubtractTowersCS::process_event(PHCompositeNode *topNode)
50 {
51  if (Verbosity() > 0)
52  std::cout << "SubtractTowersCS::process_event: entering, with _use_flow_modulation = " << _use_flow_modulation << std::endl;
53 
54  // pull out the tower containers and geometry objects at the start
55  RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
56  RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
57  RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
58 
59  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
60  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
61 
62  if (Verbosity() > 0)
63  {
64  std::cout << "SubtractTowersCS::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
65  std::cout << "SubtractTowersCS::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << std::endl;
66  std::cout << "SubtractTowersCS::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
67  }
68 
69  // these should have already been created during InitRun()
70  // note that these are the CS-versions
71  RawTowerContainer *emcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1CS");
72  RawTowerContainer *ihcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1CS");
73  RawTowerContainer *ohcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1CS");
74 
75  if (Verbosity() > 0)
76  {
77  std::cout << "SubtractTowersCS::process_event: starting with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1CS towers" << std::endl;
78  std::cout << "SubtractTowersCS::process_event: starting with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1CS towers" << std::endl;
79  std::cout << "SubtractTowersCS::process_event: starting with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1CS towers" << std::endl;
80  }
81 
82  TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub2");
83 
84  // read these in to use, even if we don't use flow modulation in the subtraction
85  float background_v2 = towerbackground->get_v2();
86  float background_Psi2 = towerbackground->get_Psi2();
87 
88  // set up constituent subtraction
89  fastjet::contrib::ConstituentSubtractor subtractor;
90 
91  // free parameter for the type of distance between particle i and
92  // ghost k. There are two options: "deltaR" or "angle" which are
93  // defined as deltaR=sqrt((y_i-y_k)^2+(phi_i-phi_k)^2) or Euclidean
94  // angle between the momenta
95  subtractor.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR);
96 
97  // free parameter for the maximal allowed distance between particle i and ghost k
98  subtractor.set_max_distance(_DeltaRmax);
99 
100  // free parameter for the distance measure (the exponent of particle
101  // pt). The larger the parameter alpha, the more are favoured the
102  // lower pt particles in the subtraction process
103  subtractor.set_alpha(_alpha);
104 
105  // free parameter for the density of ghosts. The smaller, the better
106  // - but also the computation is slower.
107  subtractor.set_ghost_area(0.01);
108 
109  // EM layer first
110  std::vector<fastjet::PseudoJet> backgroundProxies_EM;
111  std::vector<fastjet::PseudoJet> fullEvent_EM;
112  std::vector<fastjet::PseudoJet> *backgroundProxies_EM_remaining = new std::vector<fastjet::PseudoJet>;
113 
114  // create PseudoJet version of EM towers in unsubtracted event
115  RawTowerContainer::ConstRange begin_end_EM = towersEM3->getTowers();
116 
117  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
118  {
119  RawTower *tower = rtiter->second;
120 
121  double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
122  double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
123  double this_E = tower->get_energy();
124 
125  double this_pz = this_E * tanh(this_eta);
126  double this_pt = sqrt(pow(this_E, 2) - pow(this_pz, 2));
127  double this_px = this_pt * cos(this_phi);
128  double this_py = this_pt * sin(this_phi);
129 
130  fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
131 
132  fullEvent_EM.push_back(this_pj);
133  }
134 
135  // create all new towers
136  for (int eta = 0; eta < geomIH->get_etabins(); eta++)
137  {
138  for (int phi = 0; phi < geomIH->get_phibins(); phi++)
139  {
140  RawTower *new_tower = new RawTowerv1();
141  new_tower->set_energy(0);
142  emcal_towers->AddTower(eta, phi, new_tower);
143  }
144  }
145 
146  // create PseudoJet version of estimated background in EM layer
147  for (RawTowerContainer::ConstIterator rtiter = emcal_towers->getTowers().first; rtiter != emcal_towers->getTowers().second; ++rtiter)
148  {
149  RawTower *tower = rtiter->second;
150 
151  double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
152  double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
153 
154  double UE = towerbackground->get_UE(0).at(tower->get_bineta());
155  if (_use_flow_modulation)
156  {
157  UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
158  }
159 
160  double this_pz = UE * tanh(this_eta);
161  double this_pt = sqrt(pow(UE, 2) - pow(this_pz, 2));
162  double this_px = this_pt * cos(this_phi);
163  double this_py = this_pt * sin(this_phi);
164 
165  fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
166 
167  backgroundProxies_EM.push_back(this_pj);
168 
169  if (Verbosity() > 5)
170  std::cout << " SubtractTowersCS::process_event : background tower EM estimate for eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", UE = " << UE << std::endl;
171  }
172 
173  // constituent subtraction
174  std::vector<fastjet::PseudoJet> correctedEvent_EM = subtractor.do_subtraction(fullEvent_EM, backgroundProxies_EM, backgroundProxies_EM_remaining);
175 
176  if (Verbosity() > 0)
177  {
178  std::cout << " SubtractTowersCS::process_event : vector lengths fullEvent_EM = " << fullEvent_EM.size() << " , backgroundProxies_EM = " << backgroundProxies_EM.size() << " , correctedEvent_EM = " << correctedEvent_EM.size() << ", backgroundProxies_EM_remaining = " << backgroundProxies_EM_remaining->size() << std::endl;
179 
180  double E_0 = 0, E_1 = 0, E_2 = 0, E_3 = 0;
181 
182  for (unsigned int n = 0; n < fullEvent_EM.size(); n++) E_0 += fullEvent_EM.at(n).E();
183  for (unsigned int n = 0; n < backgroundProxies_EM.size(); n++) E_1 += backgroundProxies_EM.at(n).E();
184  for (unsigned int n = 0; n < correctedEvent_EM.size(); n++) E_2 += correctedEvent_EM.at(n).E();
185  for (unsigned int n = 0; n < backgroundProxies_EM_remaining->size(); n++) E_3 += backgroundProxies_EM_remaining->at(n).E();
186 
187  std::cout << "SubtractTowersCS::process_event EM : full event E - background E = " << E_0 << " - " << E_1 << " = " << E_0 - E_1 << ", subtracted E - remaining bkg E = " << E_2 << " - " << E_3 << " = " << E_2 - E_3 << std::endl;
188  }
189 
190  // load subtracted towers into grid
191 
192  for (unsigned int n = 0; n < correctedEvent_EM.size(); n++)
193  {
194  float this_eta = correctedEvent_EM.at(n).eta();
195  float this_phi = correctedEvent_EM.at(n).phi();
196  float this_E = correctedEvent_EM.at(n).E();
197 
198  // look up tower by eta / phi...
199 
200  int this_etabin = geomIH->get_etabin(this_eta);
201  int this_phibin = geomIH->get_phibin(this_phi);
202  RawTower *tower = emcal_towers->getTower(this_etabin, this_phibin);
203  tower->set_energy(this_E);
204 
205  if (Verbosity() > 5 || this_E < 0)
206  std::cout << " SubtractTowersCS::process_event : creating subtracted EM tower for eta / phi = " << this_eta << " / " << this_phi << " ( " << tower->get_binphi() << " , " << this_phibin << " ) , sub. E = " << this_E << (this_E < 0 ? " -- WARNING: negative E" : "") << std::endl;
207  }
208 
209  // IHCal layer
210  std::vector<fastjet::PseudoJet> backgroundProxies_IH;
211  std::vector<fastjet::PseudoJet> fullEvent_IH;
212  std::vector<fastjet::PseudoJet> *backgroundProxies_IH_remaining = new std::vector<fastjet::PseudoJet>;
213 
214  // create PseudoJet version of IH towers in unsubtracted event
215  RawTowerContainer::ConstRange begin_end_IH = towersIH3->getTowers();
216 
217  for (RawTowerContainer::ConstIterator rtiter = begin_end_IH.first; rtiter != begin_end_IH.second; ++rtiter)
218  {
219  RawTower *tower = rtiter->second;
220 
221  double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
222  double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
223  double this_E = tower->get_energy();
224 
225  double this_pz = this_E * tanh(this_eta);
226  double this_pt = sqrt(pow(this_E, 2) - pow(this_pz, 2));
227  double this_px = this_pt * cos(this_phi);
228  double this_py = this_pt * sin(this_phi);
229 
230  fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
231 
232  fullEvent_IH.push_back(this_pj);
233  }
234 
235  // create all new towers
236  for (int eta = 0; eta < geomIH->get_etabins(); eta++)
237  {
238  for (int phi = 0; phi < geomIH->get_phibins(); phi++)
239  {
240  RawTower *new_tower = new RawTowerv1();
241  new_tower->set_energy(0);
242  ihcal_towers->AddTower(eta, phi, new_tower);
243  }
244  }
245 
246  // create PseudoJet version of estimated background in IH layer
247  for (RawTowerContainer::ConstIterator rtiter = ihcal_towers->getTowers().first; rtiter != ihcal_towers->getTowers().second; ++rtiter)
248  {
249  RawTower *tower = rtiter->second;
250 
251  double this_eta = geomIH->get_tower_geometry(tower->get_key())->get_eta();
252  double this_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
253 
254  double UE = towerbackground->get_UE(1).at(tower->get_bineta());
255  if (_use_flow_modulation)
256  {
257  UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
258  }
259 
260  double this_pz = UE * tanh(this_eta);
261  double this_pt = sqrt(pow(UE, 2) - pow(this_pz, 2));
262  double this_px = this_pt * cos(this_phi);
263  double this_py = this_pt * sin(this_phi);
264 
265  fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
266 
267  backgroundProxies_IH.push_back(this_pj);
268 
269  if (Verbosity() > 5)
270  std::cout << " SubtractTowersCS::process_event : background tower IH estimate for eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", UE = " << UE << std::endl;
271  }
272 
273  // constituent subtraction
274  std::vector<fastjet::PseudoJet> correctedEvent_IH = subtractor.do_subtraction(fullEvent_IH, backgroundProxies_IH, backgroundProxies_IH_remaining);
275 
276  if (Verbosity() > 0)
277  {
278  std::cout << " SubtractTowersCS::process_event : vector lengths fullEvent_IH = " << fullEvent_IH.size() << " , backgroundProxies_IH = " << backgroundProxies_IH.size() << " , correctedEvent_IH = " << correctedEvent_IH.size() << ", backgroundProxies_IH_remaining = " << backgroundProxies_IH_remaining->size() << std::endl;
279 
280  double E_0 = 0, E_1 = 0, E_2 = 0, E_3 = 0;
281 
282  for (unsigned int n = 0; n < fullEvent_IH.size(); n++) E_0 += fullEvent_IH.at(n).E();
283  for (unsigned int n = 0; n < backgroundProxies_IH.size(); n++) E_1 += backgroundProxies_IH.at(n).E();
284  for (unsigned int n = 0; n < correctedEvent_IH.size(); n++) E_2 += correctedEvent_IH.at(n).E();
285  for (unsigned int n = 0; n < backgroundProxies_IH_remaining->size(); n++) E_3 += backgroundProxies_IH_remaining->at(n).E();
286 
287  std::cout << "SubtractTowersCS::process_event IH : full event E - background E = " << E_0 << " - " << E_1 << " = " << E_0 - E_1 << ", subtracted E - remaining bkg E = " << E_2 << " - " << E_3 << " = " << E_2 - E_3 << std::endl;
288  }
289 
290  // load subtracted towers into grid
291 
292  for (unsigned int n = 0; n < correctedEvent_IH.size(); n++)
293  {
294  float this_eta = correctedEvent_IH.at(n).eta();
295  float this_phi = correctedEvent_IH.at(n).phi();
296  float this_E = correctedEvent_IH.at(n).E();
297 
298  // look up tower by eta / phi...
299 
300  int this_etabin = geomIH->get_etabin(this_eta);
301  int this_phibin = geomIH->get_phibin(this_phi);
302  RawTower *tower = ihcal_towers->getTower(this_etabin, this_phibin);
303  tower->set_energy(this_E);
304 
305  if (Verbosity() > 5 || this_E < 0)
306  std::cout << " SubtractTowersCS::process_event : creating subtracted IH tower for eta / phi = " << this_eta << " / " << this_phi << " ( " << tower->get_binphi() << " , " << this_phibin << " ) , sub. E = " << this_E << (this_E < 0 ? " -- WARNING: negative E" : "") << std::endl;
307  }
308 
309  // OHCal layer
310  std::vector<fastjet::PseudoJet> backgroundProxies_OH;
311  std::vector<fastjet::PseudoJet> fullEvent_OH;
312  std::vector<fastjet::PseudoJet> *backgroundProxies_OH_remaining = new std::vector<fastjet::PseudoJet>;
313 
314  // create PseudoJet version of OH towers in unsubtracted event
315  RawTowerContainer::ConstRange begin_end_OH = towersOH3->getTowers();
316 
317  for (RawTowerContainer::ConstIterator rtiter = begin_end_OH.first; rtiter != begin_end_OH.second; ++rtiter)
318  {
319  RawTower *tower = rtiter->second;
320 
321  double this_eta = geomOH->get_tower_geometry(tower->get_key())->get_eta();
322  double this_phi = geomOH->get_tower_geometry(tower->get_key())->get_phi();
323  double this_E = tower->get_energy();
324 
325  double this_pz = this_E * tanh(this_eta);
326  double this_pt = sqrt(pow(this_E, 2) - pow(this_pz, 2));
327  double this_px = this_pt * cos(this_phi);
328  double this_py = this_pt * sin(this_phi);
329 
330  fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
331 
332  fullEvent_OH.push_back(this_pj);
333  }
334 
335  // create all new towers
336  for (int eta = 0; eta < geomOH->get_etabins(); eta++)
337  {
338  for (int phi = 0; phi < geomOH->get_phibins(); phi++)
339  {
340  RawTower *new_tower = new RawTowerv1();
341  new_tower->set_energy(0);
342  ohcal_towers->AddTower(eta, phi, new_tower);
343  }
344  }
345 
346  // create PseudoJet version of estimated background in OH layer
347  for (RawTowerContainer::ConstIterator rtiter = ohcal_towers->getTowers().first; rtiter != ohcal_towers->getTowers().second; ++rtiter)
348  {
349  RawTower *tower = rtiter->second;
350 
351  double this_eta = geomOH->get_tower_geometry(tower->get_key())->get_eta();
352  double this_phi = geomOH->get_tower_geometry(tower->get_key())->get_phi();
353 
354  double UE = towerbackground->get_UE(2).at(tower->get_bineta());
355  if (_use_flow_modulation)
356  {
357  UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
358  }
359 
360  double this_pz = UE * tanh(this_eta);
361  double this_pt = sqrt(pow(UE, 2) - pow(this_pz, 2));
362  double this_px = this_pt * cos(this_phi);
363  double this_py = this_pt * sin(this_phi);
364 
365  fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
366 
367  backgroundProxies_OH.push_back(this_pj);
368 
369  if (Verbosity() > 5)
370  std::cout << " SubtractTowersCS::process_event : background tower OH estimate for eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", UE = " << UE << std::endl;
371  }
372 
373  // constituent subtraction
374  std::vector<fastjet::PseudoJet> correctedEvent_OH = subtractor.do_subtraction(fullEvent_OH, backgroundProxies_OH, backgroundProxies_OH_remaining);
375 
376  if (Verbosity() > 0)
377  {
378  std::cout << " SubtractTowersCS::process_event : vector lengths fullEvent_OH = " << fullEvent_OH.size() << " , backgroundProxies_OH = " << backgroundProxies_OH.size() << " , correctedEvent_OH = " << correctedEvent_OH.size() << ", backgroundProxies_OH_remaining = " << backgroundProxies_OH_remaining->size() << std::endl;
379 
380  double E_0 = 0, E_1 = 0, E_2 = 0, E_3 = 0;
381 
382  for (unsigned int n = 0; n < fullEvent_OH.size(); n++) E_0 += fullEvent_OH.at(n).E();
383  for (unsigned int n = 0; n < backgroundProxies_OH.size(); n++) E_1 += backgroundProxies_OH.at(n).E();
384  for (unsigned int n = 0; n < correctedEvent_OH.size(); n++) E_2 += correctedEvent_OH.at(n).E();
385  for (unsigned int n = 0; n < backgroundProxies_OH_remaining->size(); n++) E_3 += backgroundProxies_OH_remaining->at(n).E();
386 
387  std::cout << "SubtractTowersCS::process_event OH : full event E - background E = " << E_0 << " - " << E_1 << " = " << E_0 - E_1 << ", subtracted E - remaining bkg E = " << E_2 << " - " << E_3 << " = " << E_2 - E_3 << std::endl;
388  }
389 
390  // load subtracted towers into grid
391 
392  for (unsigned int n = 0; n < correctedEvent_OH.size(); n++)
393  {
394  float this_eta = correctedEvent_OH.at(n).eta();
395  float this_phi = correctedEvent_OH.at(n).phi();
396  float this_E = correctedEvent_OH.at(n).E();
397 
398  // look up tower by eta / phi...
399 
400  int this_etabin = geomOH->get_etabin(this_eta);
401  int this_phibin = geomOH->get_phibin(this_phi);
402  RawTower *tower = ohcal_towers->getTower(this_etabin, this_phibin);
403  tower->set_energy(this_E);
404 
405  if (Verbosity() > 5 || this_E < 0)
406  std::cout << " SubtractTowersCS::process_event : creating subtracted OH tower for eta / phi = " << this_eta << " / " << this_phi << " ( " << tower->get_binphi() << " , " << this_phibin << " ) , sub. E = " << this_E << (this_E < 0 ? " -- WARNING: negative E" : "") << std::endl;
407  }
408 
409  // cleanup
410 
411  delete backgroundProxies_EM_remaining;
412  delete backgroundProxies_IH_remaining;
413  delete backgroundProxies_OH_remaining;
414 
415  if (Verbosity() > 0)
416  {
417  std::cout << "SubtractTowersCS::process_event: ending with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1CS towers" << std::endl;
418 
419  float EM_old_E = 0;
420  float EM_new_E = 0;
421  {
422  RawTowerContainer::ConstRange begin_end_EM = towersEM3->getTowers();
423  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
424  {
425  RawTower *tower = rtiter->second;
426  EM_old_E += tower->get_energy();
427  }
428  }
429  {
430  RawTowerContainer::ConstRange begin_end_EM = emcal_towers->getTowers();
431  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
432  {
433  RawTower *tower = rtiter->second;
434  EM_new_E += tower->get_energy();
435  }
436  }
437  std::cout << "SubtractTowersCS::process_event: old / new total E in EM layer = " << EM_old_E << " / " << EM_new_E << std::endl;
438 
439  std::cout << "SubtractTowersCS::process_event: ending with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1CS towers" << std::endl;
440 
441  float IH_old_E = 0;
442  float IH_new_E = 0;
443  {
444  RawTowerContainer::ConstRange begin_end_EM = towersIH3->getTowers();
445  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
446  {
447  RawTower *tower = rtiter->second;
448  IH_old_E += tower->get_energy();
449  }
450  }
451  {
452  RawTowerContainer::ConstRange begin_end_EM = ihcal_towers->getTowers();
453  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
454  {
455  RawTower *tower = rtiter->second;
456  IH_new_E += tower->get_energy();
457  }
458  }
459  std::cout << "SubtractTowersCS::process_event: old / new total E in IH layer = " << IH_old_E << " / " << IH_new_E << std::endl;
460 
461  std::cout << "SubtractTowersCS::process_event: ending with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1CS towers" << std::endl;
462 
463  float OH_old_E = 0;
464  float OH_new_E = 0;
465  {
466  RawTowerContainer::ConstRange begin_end_EM = towersOH3->getTowers();
467  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
468  {
469  RawTower *tower = rtiter->second;
470  OH_old_E += tower->get_energy();
471  }
472  }
473  {
474  RawTowerContainer::ConstRange begin_end_EM = ohcal_towers->getTowers();
475  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
476  {
477  RawTower *tower = rtiter->second;
478  OH_new_E += tower->get_energy();
479  }
480  }
481  std::cout << "SubtractTowersCS::process_event: old / new total E in OH layer = " << OH_old_E << " / " << OH_new_E << std::endl;
482  }
483 
484  if (Verbosity() > 0) std::cout << "SubtractTowersCS::process_event: exiting" << std::endl;
485 
486  return Fun4AllReturnCodes::EVENT_OK;
487 }
488 
489 int SubtractTowersCS::CreateNode(PHCompositeNode *topNode)
490 {
491  PHNodeIterator iter(topNode);
492 
493  // Looking for the DST node
494  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
495  if (!dstNode)
496  {
497  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
498  return Fun4AllReturnCodes::ABORTRUN;
499  }
500 
501  // store the new EMCal towers
502  PHCompositeNode *emcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
503  if (!emcalNode)
504  {
505  std::cout << PHWHERE << "EMCal Node note found, doing nothing." << std::endl;
506  }
507 
508  RawTowerContainer *test_emcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1CS");
509  if (!test_emcal_tower)
510  {
511  if (Verbosity() > 0) std::cout << "SubtractTowersCS::CreateNode : creating TOWER_CALIB_CEMC_RETOWER_SUB1CS node " << std::endl;
512 
513  RawTowerContainer *emcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
514  PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_towers, "TOWER_CALIB_CEMC_RETOWER_SUB1CS", "PHObject");
515  emcalNode->addNode(emcalTowerNode);
516  }
517  else
518  {
519  std::cout << "SubtractTowersCS::CreateNode : TOWER_CALIB_CEMC_RETOWER_SUB1CS already exists! " << std::endl;
520  }
521 
522  // store the new IHCal towers
523  PHCompositeNode *ihcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALIN"));
524  if (!ihcalNode)
525  {
526  std::cout << PHWHERE << "IHCal Node note found, doing nothing." << std::endl;
527  }
528 
529  RawTowerContainer *test_ihcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1CS");
530  if (!test_ihcal_tower)
531  {
532  if (Verbosity() > 0) std::cout << "SubtractTowersCS::CreateNode : creating TOWER_CALIB_HCALIN_SUB1CS node " << std::endl;
533 
534  RawTowerContainer *ihcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALIN);
535  PHIODataNode<PHObject> *ihcalTowerNode = new PHIODataNode<PHObject>(ihcal_towers, "TOWER_CALIB_HCALIN_SUB1CS", "PHObject");
536  ihcalNode->addNode(ihcalTowerNode);
537  }
538  else
539  {
540  std::cout << "SubtractTowersCS::CreateNode : TOWER_CALIB_HCALIN_SUB1CS already exists! " << std::endl;
541  }
542 
543  // store the new OHCal towers
544  PHCompositeNode *ohcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALOUT"));
545  if (!ohcalNode)
546  {
547  std::cout << PHWHERE << "OHCal Node note found, doing nothing." << std::endl;
548  }
549 
550  RawTowerContainer *test_ohcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1CS");
551  if (!test_ohcal_tower)
552  {
553  if (Verbosity() > 0) std::cout << "SubtractTowersCS::CreateNode : creating TOWER_CALIB_HCALOUT_SUB1CS node " << std::endl;
554 
555  RawTowerContainer *ohcal_towers = new RawTowerContainer(RawTowerDefs::CalorimeterId::HCALOUT);
556  PHIODataNode<PHObject> *ohcalTowerNode = new PHIODataNode<PHObject>(ohcal_towers, "TOWER_CALIB_HCALOUT_SUB1CS", "PHObject");
557  ohcalNode->addNode(ohcalTowerNode);
558  }
559  else
560  {
561  std::cout << "SubtractTowersCS::CreateNode : TOWER_CALIB_HCALOUT_SUB1CS already exists! " << std::endl;
562  }
563 
564  return Fun4AllReturnCodes::EVENT_OK;
565 }