55 , m_reco_input(
"AntiKt_Tower_r04")
56 , m_kt_input(
"Kt_Tower_r04")
57 , m_subtracted_output(
"AntiKt_Tower_r04_sub")
58 , m_etaRange(-1.0, 1.0)
59 , m_ptRange(5.0, 100.0)
61 std::cout <<
"JetMultSub::JetMultSub(const std::string &name) Calling ctor" << std::endl;
68 std::cout <<
"JetMultSub::~JetMultSub() Calling dtor" << std::endl;
75 std::cout <<
"JetMultSub::Init(PHCompositeNode *topNode) Initializing" << std::endl;
92 std::cout <<
"JetMultSub::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
99 std::cout <<
"JetMultSub::process_event(PHCompositeNode *topNode) Missing input reco jet collection" << std::endl;
107 Jet *this_jet = iter->second;
109 float this_pt = this_jet->
get_pt();
110 float this_phi = this_jet->
get_phi();
111 float this_eta = this_jet->
get_eta();
116 int nConstituents =this_jet->
size_comp();
119 float nConstituents_minus_multiplicity_correction = nConstituents - multiplicity_correction;
120 float new_pt = this_pt - rho*nConstituents_minus_multiplicity_correction;
123 float new_et = new_pt*cosh(this_eta);
124 float new_e = sqrt(pow(new_et,2) + pow(this_jet->
get_pz(),2));
125 float new_px = new_pt*cos(this_phi);
126 float new_py = new_pt*sin(this_phi);
127 float new_pz = this_jet->
get_pz();
134 new_jet->
set_e(new_e);
137 sub_jets->insert(new_jet);
142 std::cout <<
"JetMultSub::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;
143 std::cout <<
"JetMultSub::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;
171 std::cout <<
"JetMultSub::End(PHCompositeNode *topNode) This is the End..." << std::endl;
181 std::cout <<
"JetMultSub::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
191 std::cout <<
"JetMultSub::Print(const std::string &what) const Printing info for " << what << std::endl;
192 std::cout <<
"Current status: 50% working" << std::endl;
193 std::cout <<
"Errors codes: 49512323, 304, b.ss3, sos" << std::endl;
208 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
216 std::cout <<
PHWHERE <<
"ANTIKT node not found, doing nothing." << std::endl;
231 antiktNode->
addNode(subjetNode);
235 std::cout <<
"JetMultSub::CreateNode : " <<
m_subtracted_output <<
" already exists" << std::endl;
250 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Estimating rho" << std::endl;
256 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Missing input kT jet collection" << std::endl;
263 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Found input kT jet collection" << std::endl;
264 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Number of kT jets: " << kt_jets->
size() << std::endl;
268 int hardest_kT_jet_indices[2] = {-1,-1};
269 float hardest_kT_jet_pts[2] = {-1.0,-1.0};
272 Jet *jet = iter->second;
279 if (pt > hardest_kT_jet_pts[0])
281 hardest_kT_jet_pts[1] = hardest_kT_jet_pts[0];
282 hardest_kT_jet_indices[1] = hardest_kT_jet_indices[0];
283 hardest_kT_jet_pts[0] =
pt;
284 hardest_kT_jet_indices[0] = jet->
get_id();
286 else if (pt > hardest_kT_jet_pts[1])
288 hardest_kT_jet_pts[1] =
pt;
289 hardest_kT_jet_indices[1] = jet->
get_id();
298 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Hardest kT jet indices: " << hardest_kT_jet_indices[0] <<
", " << hardest_kT_jet_indices[1] << std::endl;
299 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Hardest kT jet et: " << hardest_kT_jet_pts[0] <<
", " << hardest_kT_jet_pts[1] << std::endl;
303 std::vector<float> et_over_nConstituents;
306 Jet *jet = iter->second;
315 if (
id != hardest_kT_jet_indices[0] &&
id != hardest_kT_jet_indices[1])
317 et_over_nConstituents.push_back(pt/nConstituents);
323 float median_et_over_nConstituents = 0.0;
324 if (et_over_nConstituents.size() > 0)
326 std::sort(et_over_nConstituents.begin(), et_over_nConstituents.end());
327 median_et_over_nConstituents = et_over_nConstituents[et_over_nConstituents.size()/2];
333 std::cout <<
"JetMultSub::EstimateRho(PHCompositeNode *topNode) Median et over nConstituents: " << median_et_over_nConstituents << std::endl;
337 return median_et_over_nConstituents;
346 std::cout <<
"JetMultSub::GetMultiplicityCorrection(float pt) Getting multiplicity correction" << std::endl;
351 std::vector<float> pt_reco_bins = {0.0,100.0};
352 std::vector<float> multiplicity_corrections = {0.0};
355 int pt_size = pt_reco_bins.size();
356 for (
int i = 0;
i < pt_size;
i++){
357 if (pt > pt_reco_bins[
i]){
366 correction = multiplicity_corrections[pt_bin];