5 #include <calobase/RawTower.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerGeom.h>
8 #include <calobase/RawTowerGeomContainer.h>
35 , _use_flow_modulation(
false)
49 std::cout <<
"CopyAndSubtractJets::process_event: entering, with _use_flow_modulation = " <<
_use_flow_modulation << std::endl;
52 RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_CEMC_RETOWER");
53 RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALIN");
54 RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode,
"TOWER_CALIB_HCALOUT");
56 RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALIN");
57 RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode,
"TOWERGEOM_HCALOUT");
60 JetMap *unsub_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_HIRecoSeedsRaw_r02");
61 JetMap *sub_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_HIRecoSeedsSub_r02");
63 TowerBackground *background = findNode::getClass<TowerBackground>(topNode,
"TowerBackground_Sub1");
65 std::vector<float> background_UE_0 = background->
get_UE(0);
66 std::vector<float> background_UE_1 = background->get_UE(1);
67 std::vector<float> background_UE_2 = background->get_UE(2);
69 float background_v2 = background->get_v2();
70 float background_Psi2 = background->get_Psi2();
74 std::cout <<
"CopyAndSubtractJets::process_event: entering with # unsubtracted jets = " << unsub_jets->size() << std::endl;
75 std::cout <<
"CopyAndSubtractJets::process_event: entering with # subtracted jets = " << sub_jets->size() << std::endl;
80 for (
JetMap::Iter iter = unsub_jets->begin(); iter != unsub_jets->end(); ++iter)
82 Jet *this_jet = iter->second;
84 float this_pt = this_jet->
get_pt();
85 float this_phi = this_jet->
get_phi();
86 float this_eta = this_jet->
get_eta();
90 float new_total_px = 0;
91 float new_total_py = 0;
92 float new_total_pz = 0;
93 float new_total_e = 0;
98 std::cout <<
"CopyAndSubtractJets::process_event: unsubtracted jet with pt / eta / phi = " << this_pt <<
" / " << this_eta <<
" / " << this_phi << std::endl;
111 double comp_background = 0;
113 if ((*comp).first == 5)
115 tower = towersIH3->getTower((*comp).second);
116 tower_geom = geomIH->get_tower_geometry(tower->
get_key());
118 comp_ieta = geomIH->get_etabin(tower_geom->
get_eta());
119 comp_background = background_UE_1.at(comp_ieta);
121 else if ((*comp).first == 7)
123 tower = towersOH3->getTower((*comp).second);
124 tower_geom = geomOH->get_tower_geometry(tower->
get_key());
126 comp_ieta = geomOH->get_etabin(tower_geom->
get_eta());
127 comp_background = background_UE_2.at(comp_ieta);
129 else if ((*comp).first == 13)
131 tower = towersEM3->
getTower((*comp).second);
132 tower_geom = geomIH->get_tower_geometry(tower->
get_key());
134 comp_ieta = geomIH->get_etabin(tower_geom->
get_eta());
135 comp_background = background_UE_0.at(comp_ieta);
142 comp_eta = tower_geom->
get_eta();
143 comp_phi = tower_geom->
get_phi();
148 std::cout <<
"CopyAndSubtractJets::process_event: --> constituent in layer " << (*comp).first <<
", has unsub E = " << comp_e <<
", is at ieta #" << comp_ieta <<
", and has UE = " << comp_background << std::endl;
154 comp_background = comp_background * (1 + 2 * background_v2 * cos(2 * (comp_phi - background_Psi2)));
156 std::cout <<
"CopyAndSubtractJets::process_event: --> --> flow mod, at phi = " << comp_phi <<
", v2 and Psi2 are = " << background_v2 <<
" , " << background_Psi2 <<
", UE after modulation = " << comp_background << std::endl;
160 double comp_sub_e = comp_e - comp_background;
164 double comp_px = comp_sub_e / cosh(comp_eta) * cos(comp_phi);
165 double comp_py = comp_sub_e / cosh(comp_eta) * sin(comp_phi);
166 double comp_pz = comp_sub_e * tanh(comp_eta);
168 new_total_px += comp_px;
169 new_total_py += comp_py;
170 new_total_pz += comp_pz;
171 new_total_e += comp_sub_e;
174 new_jet->
set_px(new_total_px);
175 new_jet->
set_py(new_total_py);
176 new_jet->
set_pz(new_total_pz);
177 new_jet->
set_e(new_total_e);
180 sub_jets->insert(new_jet);
184 std::cout <<
"CopyAndSubtractJets::process_event: old jet #" << ijet <<
", old px / py / pz / e = " << this_jet->
get_px() <<
" / " << this_jet->
get_py() <<
" / " << this_jet->
get_pz() <<
" / " << this_jet->
get_e() << std::endl;
185 std::cout <<
"CopyAndSubtractJets::process_event: new jet #" << ijet <<
", new px / py / pz / e = " << new_jet->
get_px() <<
" / " << new_jet->
get_py() <<
" / " << new_jet->
get_pz() <<
" / " << new_jet->
get_e() << std::endl;
193 std::cout <<
"CopyAndSubtractJets::process_event: exiting with # subtracted jets = " << sub_jets->size() << std::endl;
212 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
220 std::cout <<
PHWHERE <<
"ANTIKT node not found, doing nothing." << std::endl;
227 std::cout <<
PHWHERE <<
"TOWER node not found, doing nothing." << std::endl;
231 JetMap *test_jets = findNode::getClass<JetMap>(topNode,
"AntiKt_Tower_HIRecoSeedsSub_r02");
234 if (
Verbosity() > 0) std::cout <<
"CopyAndSubtractJets::CreateNode : creating AntiKt_Tower_HIRecoSeedsSub_r02 node " << std::endl;
238 towerNode->
addNode(subjetNode);
242 std::cout <<
"CopyAndSubtractJets::CreateNode : AntiKt_Tower_HIRecoSeedsSub_r02 already exists! " << std::endl;