ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SubtractTowers.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SubtractTowers.cc
1 #include "SubtractTowers.h"
2 
3 #include "TowerBackground.h"
4 
5 // sPHENIX includes
6 #include <calobase/RawTower.h>
7 #include <calobase/RawTowerContainer.h>
8 #include <calobase/RawTowerDefs.h>
9 #include <calobase/RawTowerGeom.h>
10 #include <calobase/RawTowerGeomContainer.h>
11 #include <calobase/RawTowerv1.h>
12 
14 #include <fun4all/SubsysReco.h>
15 
16 #include <phool/PHCompositeNode.h>
17 #include <phool/PHIODataNode.h>
18 #include <phool/PHNode.h>
19 #include <phool/PHNodeIterator.h>
20 #include <phool/PHObject.h>
21 #include <phool/getClass.h>
22 #include <phool/phool.h>
23 
24 // standard includes
25 #include <cmath>
26 #include <iostream>
27 #include <map>
28 #include <utility>
29 #include <vector>
30 
32  : SubsysReco(name)
33  , _use_flow_modulation(false)
34 {
35 }
36 
38 {
39  CreateNode(topNode);
40 
42 }
43 
45 {
46  if (Verbosity() > 0)
47  std::cout << "SubtractTowers::process_event: entering, with _use_flow_modulation = " << _use_flow_modulation << std::endl;
48 
49  // pull out the tower containers and geometry objects at the start
50  RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
51  RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
52  RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
53 
54  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
55  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
56 
57  if (Verbosity() > 0)
58  {
59  std::cout << "SubtractTowers::process_event: " << towersEM3->size() << " TOWER_CALIB_CEMC_RETOWER towers" << std::endl;
60  std::cout << "SubtractTowers::process_event: " << towersIH3->size() << " TOWER_CALIB_HCALIN towers" << std::endl;
61  std::cout << "SubtractTowers::process_event: " << towersOH3->size() << " TOWER_CALIB_HCALOUT towers" << std::endl;
62  }
63 
64  // these should have already been created during InitRun()
65  RawTowerContainer *emcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
66  RawTowerContainer *ihcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
67  RawTowerContainer *ohcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
68 
69  if (Verbosity() > 0)
70  {
71  std::cout << "SubtractTowers::process_event: starting with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1 towers" << std::endl;
72  std::cout << "SubtractTowers::process_event: starting with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1 towers" << std::endl;
73  std::cout << "SubtractTowers::process_event: starting with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1 towers" << std::endl;
74  }
75 
76  TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub2");
77 
78  // read these in to use, even if we don't use flow modulation in the subtraction
79  float background_v2 = towerbackground->get_v2();
80  float background_Psi2 = towerbackground->get_Psi2();
81 
82  // EMCal
83 
84  // replicate existing towers
85  RawTowerContainer::ConstRange begin_end_EM = towersEM3->getTowers();
86  for (RawTowerContainer::ConstIterator rtiter = begin_end_EM.first; rtiter != begin_end_EM.second; ++rtiter)
87  {
88  RawTower *tower = rtiter->second;
89 
90  int this_etabin = tower->get_bineta();
91  int this_phibin = tower->get_binphi();
92  float this_E = tower->get_energy();
93 
94  RawTower *new_tower = new RawTowerv1();
95  new_tower->set_energy(this_E);
96  emcal_towers->AddTower(this_etabin, this_phibin, new_tower);
97  }
98 
99  // now fill in additional towers with zero energy to fill out the full grid
100  // but note: after retowering, all of these should already exist...
101  for (int eta = 0; eta < geomIH->get_etabins(); eta++)
102  {
103  for (int phi = 0; phi < geomIH->get_phibins(); phi++)
104  {
105  if (!emcal_towers->getTower(eta, phi))
106  {
107  RawTower *new_tower = new RawTowerv1();
108  new_tower->set_energy(0);
109  emcal_towers->AddTower(eta, phi, new_tower);
110  }
111  }
112  }
113 
114  // update towers for background subtraction...
115  for (RawTowerContainer::ConstIterator rtiter = emcal_towers->getTowers().first; rtiter != emcal_towers->getTowers().second; ++rtiter)
116  {
117  RawTower *tower = rtiter->second;
118  float raw_energy = tower->get_energy();
119  float UE = towerbackground->get_UE(0).at(tower->get_bineta());
121  {
122  float tower_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
123  UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
124  }
125  float new_energy = raw_energy - UE;
126  tower->set_energy(new_energy);
127  if (Verbosity() > 5)
128  std::cout << " SubtractTowers::process_event : EMCal tower at eta / phi = " << tower->get_bineta() << " / " << tower->get_binphi() << ", pre-sub / after-sub E = " << raw_energy << " / " << tower->get_energy() << std::endl;
129  }
130 
131  // IHCal
132 
133  // replicate existing towers
134  RawTowerContainer::ConstRange begin_end_IH = towersIH3->getTowers();
135  for (RawTowerContainer::ConstIterator rtiter = begin_end_IH.first; rtiter != begin_end_IH.second; ++rtiter)
136  {
137  RawTower *tower = rtiter->second;
138  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
139 
140  int this_etabin = geomIH->get_etabin(tower_geom->get_eta());
141  int this_phibin = geomIH->get_phibin(tower_geom->get_phi());
142  float this_E = tower->get_energy();
143 
144  RawTower *new_tower = new RawTowerv1();
145  new_tower->set_energy(this_E);
146  ihcal_towers->AddTower(this_etabin, this_phibin, new_tower);
147  }
148 
149  // now fill in additional towers with zero energy to fill out the full grid
150  for (int eta = 0; eta < geomIH->get_etabins(); eta++)
151  {
152  for (int phi = 0; phi < geomIH->get_phibins(); phi++)
153  {
154  if (!ihcal_towers->getTower(eta, phi))
155  {
156  RawTower *new_tower = new RawTowerv1();
157  new_tower->set_energy(0);
158  ihcal_towers->AddTower(eta, phi, new_tower);
159  }
160  }
161  }
162 
163  // update towers for background subtraction...
164  for (RawTowerContainer::ConstIterator rtiter = ihcal_towers->getTowers().first; rtiter != ihcal_towers->getTowers().second; ++rtiter)
165  {
166  RawTower *tower = rtiter->second;
167  float raw_energy = tower->get_energy();
168  float UE = towerbackground->get_UE(1).at(tower->get_bineta());
170  {
171  float tower_phi = geomIH->get_tower_geometry(tower->get_key())->get_phi();
172  UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
173  }
174  float new_energy = raw_energy - UE;
175  tower->set_energy(new_energy);
176  }
177 
178  // OHCal
179 
180  // replicate existing towers
181  RawTowerContainer::ConstRange begin_end_OH = towersOH3->getTowers();
182  for (RawTowerContainer::ConstIterator rtiter = begin_end_OH.first; rtiter != begin_end_OH.second; ++rtiter)
183  {
184  RawTower *tower = rtiter->second;
185  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
186 
187  int this_etabin = geomOH->get_etabin(tower_geom->get_eta());
188  int this_phibin = geomOH->get_phibin(tower_geom->get_phi());
189  float this_E = tower->get_energy();
190 
191  RawTower *new_tower = new RawTowerv1();
192  new_tower->set_energy(this_E);
193  ohcal_towers->AddTower(this_etabin, this_phibin, new_tower);
194  }
195 
196  // now fill in additional towers with zero energy to fill out the full grid
197  for (int eta = 0; eta < geomOH->get_etabins(); eta++)
198  {
199  for (int phi = 0; phi < geomOH->get_phibins(); phi++)
200  {
201  if (!ohcal_towers->getTower(eta, phi))
202  {
203  RawTower *new_tower = new RawTowerv1();
204  new_tower->set_energy(0);
205  ohcal_towers->AddTower(eta, phi, new_tower);
206  }
207  }
208  }
209 
210  // update towers for background subtraction...
211  for (RawTowerContainer::ConstIterator rtiter = ohcal_towers->getTowers().first; rtiter != ohcal_towers->getTowers().second; ++rtiter)
212  {
213  RawTower *tower = rtiter->second;
214  float raw_energy = tower->get_energy();
215  float UE = towerbackground->get_UE(2).at(tower->get_bineta());
217  {
218  float tower_phi = geomOH->get_tower_geometry(tower->get_key())->get_phi();
219  UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
220  }
221  float new_energy = raw_energy - UE;
222  tower->set_energy(new_energy);
223  }
224 
225  if (Verbosity() > 0)
226  {
227  std::cout << "SubtractTowers::process_event: ending with " << emcal_towers->size() << " TOWER_CALIB_CEMC_RETOWER_SUB1 towers" << std::endl;
228  std::cout << "SubtractTowers::process_event: ending with " << ihcal_towers->size() << " TOWER_CALIB_HCALIN_SUB1 towers" << std::endl;
229  std::cout << "SubtractTowers::process_event: ending with " << ohcal_towers->size() << " TOWER_CALIB_HCALOUT_SUB1 towers" << std::endl;
230  }
231 
232  if (Verbosity() > 0) std::cout << "SubtractTowers::process_event: exiting" << std::endl;
233 
235 }
236 
238 {
239  PHNodeIterator iter(topNode);
240 
241  // Looking for the DST node
242  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
243  if (!dstNode)
244  {
245  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
247  }
248 
249  // store the new EMCal towers
250  PHCompositeNode *emcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
251  if (!emcalNode)
252  {
253  std::cout << PHWHERE << "EMCal Node note found, doing nothing." << std::endl;
254  }
255 
256  RawTowerContainer *test_emcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
257  if (!test_emcal_tower)
258  {
259  if (Verbosity() > 0) std::cout << "SubtractTowers::CreateNode : creating TOWER_CALIB_CEMC_RETOWER_SUB1 node " << std::endl;
260 
262  PHIODataNode<PHObject> *emcalTowerNode = new PHIODataNode<PHObject>(emcal_towers, "TOWER_CALIB_CEMC_RETOWER_SUB1", "PHObject");
263  emcalNode->addNode(emcalTowerNode);
264  }
265  else
266  {
267  std::cout << "SubtractTowers::CreateNode : TOWER_CALIB_CEMC_RETOWER_SUB1 already exists! " << std::endl;
268  }
269 
270  // store the new IHCal towers
271  PHCompositeNode *ihcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALIN"));
272  if (!ihcalNode)
273  {
274  std::cout << PHWHERE << "IHCal Node note found, doing nothing." << std::endl;
275  }
276 
277  RawTowerContainer *test_ihcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
278  if (!test_ihcal_tower)
279  {
280  if (Verbosity() > 0) std::cout << "SubtractTowers::CreateNode : creating TOWER_CALIB_HCALIN_SUB1 node " << std::endl;
281 
283  PHIODataNode<PHObject> *ihcalTowerNode = new PHIODataNode<PHObject>(ihcal_towers, "TOWER_CALIB_HCALIN_SUB1", "PHObject");
284  ihcalNode->addNode(ihcalTowerNode);
285  }
286  else
287  {
288  std::cout << "SubtractTowers::CreateNode : TOWER_CALIB_HCALIN_SUB1 already exists! " << std::endl;
289  }
290 
291  // store the new OHCal towers
292  PHCompositeNode *ohcalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "HCALOUT"));
293  if (!ohcalNode)
294  {
295  std::cout << PHWHERE << "OHCal Node note found, doing nothing." << std::endl;
296  }
297 
298  RawTowerContainer *test_ohcal_tower = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
299  if (!test_ohcal_tower)
300  {
301  if (Verbosity() > 0) std::cout << "SubtractTowers::CreateNode : creating TOWER_CALIB_HCALOUT_SUB1 node " << std::endl;
302 
304  PHIODataNode<PHObject> *ohcalTowerNode = new PHIODataNode<PHObject>(ohcal_towers, "TOWER_CALIB_HCALOUT_SUB1", "PHObject");
305  ohcalNode->addNode(ohcalTowerNode);
306  }
307  else
308  {
309  std::cout << "SubtractTowers::CreateNode : TOWER_CALIB_HCALOUT_SUB1 already exists! " << std::endl;
310  }
311 
313 }