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>
23 #include <fastjet/PseudoJet.hh>
24 #include <fastjet/contrib/ConstituentSubtractor.hh>
36 , _use_flow_modulation(
false)
52 std::cout <<
"SubtractTowersCS::process_event: entering, with _use_flow_modulation = " <<
_use_flow_modulation << std::endl;
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");
59 RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALIN");
60 RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALOUT");
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;
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");
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;
82 TowerBackground *towerbackground = findNode::getClass<TowerBackground>(topNode,
"TowerBackground_Sub2");
85 float background_v2 = towerbackground->
get_v2();
86 float background_Psi2 = towerbackground->
get_Psi2();
89 fastjet::contrib::ConstituentSubtractor subtractor;
103 subtractor.set_alpha(
_alpha);
107 subtractor.set_ghost_area(0.01);
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>;
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();
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);
130 fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
132 fullEvent_EM.push_back(this_pj);
136 for (
int eta = 0;
eta < geomIH->get_etabins();
eta++)
138 for (
int phi = 0;
phi < geomIH->get_phibins();
phi++)
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();
157 UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
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);
165 fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
167 backgroundProxies_EM.push_back(this_pj);
170 std::cout <<
" SubtractTowersCS::process_event : background tower EM estimate for eta / phi = " << tower->
get_bineta() <<
" / " << tower->
get_binphi() <<
", UE = " << UE << std::endl;
174 std::vector<fastjet::PseudoJet> correctedEvent_EM = subtractor.do_subtraction(fullEvent_EM, backgroundProxies_EM, backgroundProxies_EM_remaining);
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;
180 double E_0 = 0, E_1 = 0, E_2 = 0, E_3 = 0;
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();
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;
192 for (
unsigned int n = 0;
n < correctedEvent_EM.size();
n++)
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();
200 int this_etabin = geomIH->get_etabin(this_eta);
201 int this_phibin = geomIH->get_phibin(this_phi);
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;
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>;
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();
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);
230 fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
232 fullEvent_IH.push_back(this_pj);
236 for (
int eta = 0;
eta < geomIH->get_etabins();
eta++)
238 for (
int phi = 0;
phi < geomIH->get_phibins();
phi++)
242 ihcal_towers->AddTower(
eta,
phi, new_tower);
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();
257 UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
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);
265 fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
267 backgroundProxies_IH.push_back(this_pj);
270 std::cout <<
" SubtractTowersCS::process_event : background tower IH estimate for eta / phi = " << tower->
get_bineta() <<
" / " << tower->
get_binphi() <<
", UE = " << UE << std::endl;
274 std::vector<fastjet::PseudoJet> correctedEvent_IH = subtractor.do_subtraction(fullEvent_IH, backgroundProxies_IH, backgroundProxies_IH_remaining);
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;
280 double E_0 = 0, E_1 = 0, E_2 = 0, E_3 = 0;
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();
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;
292 for (
unsigned int n = 0;
n < correctedEvent_IH.size();
n++)
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();
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);
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;
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>;
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();
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);
330 fastjet::PseudoJet this_pj(this_px, this_py, this_pz, this_E);
332 fullEvent_OH.push_back(this_pj);
336 for (
int eta = 0;
eta < geomOH->get_etabins();
eta++)
338 for (
int phi = 0;
phi < geomOH->get_phibins();
phi++)
342 ohcal_towers->AddTower(
eta,
phi, new_tower);
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();
357 UE = UE * (1 + 2 * background_v2 * cos(2 * (this_phi - background_Psi2)));
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);
365 fastjet::PseudoJet this_pj(this_px, this_py, this_pz, UE);
367 backgroundProxies_OH.push_back(this_pj);
370 std::cout <<
" SubtractTowersCS::process_event : background tower OH estimate for eta / phi = " << tower->
get_bineta() <<
" / " << tower->
get_binphi() <<
", UE = " << UE << std::endl;
374 std::vector<fastjet::PseudoJet> correctedEvent_OH = subtractor.do_subtraction(fullEvent_OH, backgroundProxies_OH, backgroundProxies_OH_remaining);
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;
380 double E_0 = 0, E_1 = 0, E_2 = 0, E_3 = 0;
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();
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;
392 for (
unsigned int n = 0;
n < correctedEvent_OH.size();
n++)
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();
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);
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;
411 delete backgroundProxies_EM_remaining;
412 delete backgroundProxies_IH_remaining;
413 delete backgroundProxies_OH_remaining;
417 std::cout <<
"SubtractTowersCS::process_event: ending with " << emcal_towers->
size() <<
" TOWER_CALIB_CEMC_RETOWER_SUB1CS towers" << std::endl;
437 std::cout <<
"SubtractTowersCS::process_event: old / new total E in EM layer = " << EM_old_E <<
" / " << EM_new_E << std::endl;
439 std::cout <<
"SubtractTowersCS::process_event: ending with " << ihcal_towers->size() <<
" TOWER_CALIB_HCALIN_SUB1CS towers" << std::endl;
459 std::cout <<
"SubtractTowersCS::process_event: old / new total E in IH layer = " << IH_old_E <<
" / " << IH_new_E << std::endl;
461 std::cout <<
"SubtractTowersCS::process_event: ending with " << ohcal_towers->size() <<
" TOWER_CALIB_HCALOUT_SUB1CS towers" << std::endl;
481 std::cout <<
"SubtractTowersCS::process_event: old / new total E in OH layer = " << OH_old_E <<
" / " << OH_new_E << std::endl;
484 if (
Verbosity() > 0) std::cout <<
"SubtractTowersCS::process_event: exiting" << std::endl;
497 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
505 std::cout <<
PHWHERE <<
"EMCal Node note found, doing nothing." << std::endl;
508 RawTowerContainer *test_emcal_tower = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_CEMC_RETOWER_SUB1CS");
509 if (!test_emcal_tower)
511 if (
Verbosity() > 0) std::cout <<
"SubtractTowersCS::CreateNode : creating TOWER_CALIB_CEMC_RETOWER_SUB1CS node " << std::endl;
515 emcalNode->
addNode(emcalTowerNode);
519 std::cout <<
"SubtractTowersCS::CreateNode : TOWER_CALIB_CEMC_RETOWER_SUB1CS already exists! " << std::endl;
526 std::cout <<
PHWHERE <<
"IHCal Node note found, doing nothing." << std::endl;
529 RawTowerContainer *test_ihcal_tower = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALIN_SUB1CS");
530 if (!test_ihcal_tower)
532 if (
Verbosity() > 0) std::cout <<
"SubtractTowersCS::CreateNode : creating TOWER_CALIB_HCALIN_SUB1CS node " << std::endl;
536 ihcalNode->
addNode(ihcalTowerNode);
540 std::cout <<
"SubtractTowersCS::CreateNode : TOWER_CALIB_HCALIN_SUB1CS already exists! " << std::endl;
547 std::cout <<
PHWHERE <<
"OHCal Node note found, doing nothing." << std::endl;
550 RawTowerContainer *test_ohcal_tower = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALOUT_SUB1CS");
551 if (!test_ohcal_tower)
553 if (
Verbosity() > 0) std::cout <<
"SubtractTowersCS::CreateNode : creating TOWER_CALIB_HCALOUT_SUB1CS node " << std::endl;
557 ohcalNode->
addNode(ohcalTowerNode);
561 std::cout <<
"SubtractTowersCS::CreateNode : TOWER_CALIB_HCALOUT_SUB1CS already exists! " << std::endl;