Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JetMultSub.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetMultSub.cc
1 /*
2 JetMultSub SubsysReco module.
3 
4  Description:
5  This is a module that takes in reconstructed jets and subtracts the background using the multiplicity method.
6  Estimates rho using the median of the average pt in a kT jet with radius 0.4.
7  Corrects the multiplicity using the average truth multiplicity of a a jet with a given pt_reco.
8  Subtracts the background from the jet pt using the formula:
9  pt_sub = pt_reco - rho*(N_reco - <N_signal>(pt_reco))
10 
11  Author: Tanner Mengel
12 
13  Date: 10/04/2023
14 
15 */
16 
17 //____________________________________________________________________________..
18 // JetMultSub.h
19 #include "JetMultSub.h"
20 //____________________________________________________________________________..
21 // Fun4All includes
23 #include <fun4all/SubsysReco.h>
24 //____________________________________________________________________________..
25 // PHOOL includes
26 #include <phool/PHCompositeNode.h>
27 #include <phool/getClass.h>
28 #include <phool/PHCompositeNode.h>
29 #include <phool/PHIODataNode.h>
30 #include <phool/PHNode.h>
31 #include <phool/PHNodeIterator.h>
32 #include <phool/PHObject.h>
33 #include <phool/phool.h>
34 //____________________________________________________________________________..
35 // Jet includes
36 #include <jetbase/Jet.h>
37 #include <jetbase/JetMap.h>
38 #include <jetbase/JetMapv1.h>
39 #include <jetbase/Jetv1.h>
40 //____________________________________________________________________________..
41 // standard includes
42 #include <cmath>
43 #include <iostream>
44 #include <map>
45 #include <vector>
46 #include <utility>
47 //____________________________________________________________________________..
48 // Starting function definitions...
49 
50 
51 //____________________________________________________________________________..
52 // constructor. Sets default values
54  SubsysReco(name)
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)
60 {
61  std::cout << "JetMultSub::JetMultSub(const std::string &name) Calling ctor" << std::endl;
62 }
63 
64 //____________________________________________________________________________..
65 // destructor. Does nothing
67 {
68  std::cout << "JetMultSub::~JetMultSub() Calling dtor" << std::endl;
69 }
70 
71 //____________________________________________________________________________..
72 // init. Does nothing
74 {
75  std::cout << "JetMultSub::Init(PHCompositeNode *topNode) Initializing" << std::endl;
77 }
78 
79 //____________________________________________________________________________..
80 // init run. Very important! Creates output node if it doesn't exist
82 {
83  CreateNode(topNode);
85 }
86 
87 //____________________________________________________________________________..
88 // process event. This is where the magic happens
90 {
91  if (Verbosity() > 0){
92  std::cout << "JetMultSub::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
93  }
94 
95  // get unsubtracted and subtracted jet collections
96  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, m_reco_input);
97  JetMap *sub_jets = findNode::getClass<JetMap>(topNode, m_subtracted_output);
98  if(!reco_jets){
99  std::cout << "JetMultSub::process_event(PHCompositeNode *topNode) Missing input reco jet collection" << std::endl;
101  }
102 
103  // loop over reco jets and subtract background
104  int ijet = 0;
105  for (JetMap::Iter iter = reco_jets->begin(); iter != reco_jets->end(); ++iter)
106  {
107  Jet *this_jet = iter->second;
108 
109  float this_pt = this_jet->get_pt();
110  float this_phi = this_jet->get_phi();
111  float this_eta = this_jet->get_eta();
112 
113  Jet *new_jet = new Jetv1();
114 
115  // subtract background
116  int nConstituents =this_jet->size_comp();
117  float rho = EstimateRho(topNode);
118  float multiplicity_correction = GetMultiplicityCorrection(this_pt);
119  float nConstituents_minus_multiplicity_correction = nConstituents - multiplicity_correction;
120  float new_pt = this_pt - rho*nConstituents_minus_multiplicity_correction;
121 
122  // calculate new px, py, pz, e
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();
128 
129 
130  // set new jet properties
131  new_jet->set_px(new_px);
132  new_jet->set_py(new_py);
133  new_jet->set_pz(new_pz);
134  new_jet->set_e(new_e);
135 
136  // add new jet to subtracted jet collection
137  sub_jets->insert(new_jet);
138 
139  // debug info
140  if (Verbosity() > 1 && this_pt > 5)
141  {
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;
144  }
145 
146  ijet++;
147 
148  }
149 
151 }
152 
153 //____________________________________________________________________________..
154 // reset event. Does nothing
156 {
158 }
159 
160 //____________________________________________________________________________..
161 // end run. You guessed it! does nothing
163 {
165 }
166 
167 //____________________________________________________________________________..'
168 // this is the end. Does nothing
170 {
171  std::cout << "JetMultSub::End(PHCompositeNode *topNode) This is the End..." << std::endl;
173 }
174 
175 //____________________________________________________________________________..
176 // reset. Does nothing
178 {
179  if (Verbosity() > 0)
180  {
181  std::cout << "JetMultSub::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
182  }
184 }
185 
186 //____________________________________________________________________________..
187 // print. Prints meaningless message to screen
188 void JetMultSub::Print(const std::string &what) const
189 {
190  if (Verbosity() > 0){
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;
194  }
195 }
196 
197 //____________________________________________________________________________..
198 // create node. This is actually important. Creates the node if it doesn't exist
200 {
201  // iterate to the DST node and add the jet collection if necessary
202  PHNodeIterator iter(topNode);
203 
204  // Looking for the DST node
205  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
206  if (!dstNode)
207  {
208  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
210  }
211 
212  // Looking for the ANTIKT node
213  PHCompositeNode *antiktNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "ANTIKT"));
214  if (!antiktNode)
215  {
216  std::cout << PHWHERE << "ANTIKT node not found, doing nothing." << std::endl;
217  }
218 
219  // store the new jet collection if it doesn't exist already
220  JetMap *test_jets;
221  test_jets = findNode::getClass<JetMap>(topNode, m_subtracted_output);
222  if(!test_jets)
223  {
224  JetMap *sub_jets = new JetMapv1();
225  PHIODataNode<PHObject> *subjetNode ;
226  if (Verbosity() > 0)
227  {
228  std::cout << "JetMultSub::CreateNode : creating " << m_subtracted_output << std::endl;
229  }
230  subjetNode = new PHIODataNode<PHObject>(sub_jets, m_subtracted_output, "PHObject");
231  antiktNode->addNode(subjetNode);
232  }
233  else
234  {
235  std::cout << "JetMultSub::CreateNode : " << m_subtracted_output << " already exists" << std::endl;
236  }
237 
238  // all done!
239 
241 
242 }
243 
244 //____________________________________________________________________________..
245 // estimate rho. Gets the event rho value from r=0.4 kT jets
247 {
248 
249  if (Verbosity() > 0){
250  std::cout << "JetMultSub::EstimateRho(PHCompositeNode *topNode) Estimating rho" << std::endl;
251  }
252 
253  // get the kT jet collection
254  JetMap *kt_jets = findNode::getClass<JetMap>(topNode, m_kt_input);
255  if(!kt_jets){
256  std::cout << "JetMultSub::EstimateRho(PHCompositeNode *topNode) Missing input kT jet collection" << std::endl;
257  return -1.0;
258  }
259 
260  // debug info
261  if (Verbosity() > 0)
262  {
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;
265  }
266 
267  // loop over kT jets and find the 2 hardest jets in the eta range
268  int hardest_kT_jet_indices[2] = {-1,-1};
269  float hardest_kT_jet_pts[2] = {-1.0,-1.0};
270  for (JetMap::Iter iter = kt_jets->begin(); iter != kt_jets->end(); ++iter)
271  {
272  Jet *jet = iter->second;
273  float pt = jet->get_pt();
274  float eta = jet->get_eta();
275  // float et = jet->get_et();
276 
277  if (eta > m_etaRange.first && eta < m_etaRange.second)
278  {
279  if (pt > hardest_kT_jet_pts[0])
280  {
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();
285  }
286  else if (pt > hardest_kT_jet_pts[1])
287  {
288  hardest_kT_jet_pts[1] = pt;
289  hardest_kT_jet_indices[1] = jet->get_id();
290  }
291  }
292  }
293 
294 
295  // debug info
296  if(Verbosity() > 0)
297  {
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;
300  }
301 
302  // get median value of the average et of constituents in all but the 2 hardest jets
303  std::vector<float> et_over_nConstituents;
304  for (JetMap::Iter iter = kt_jets->begin(); iter != kt_jets->end(); ++iter)
305  {
306  Jet *jet = iter->second;
307  float pt = jet->get_pt();
308  float eta = jet->get_eta();
309  // float et = jet->get_et();
310  int nConstituents = jet->size_comp();
311 
312  if (eta > m_etaRange.first && eta < m_etaRange.second)
313  {
314  int id = jet->get_id();
315  if (id != hardest_kT_jet_indices[0] && id != hardest_kT_jet_indices[1])
316  {
317  et_over_nConstituents.push_back(pt/nConstituents);
318  }
319  }
320  }
321 
322  // get median value
323  float median_et_over_nConstituents = 0.0;
324  if (et_over_nConstituents.size() > 0)
325  {
326  std::sort(et_over_nConstituents.begin(), et_over_nConstituents.end());
327  median_et_over_nConstituents = et_over_nConstituents[et_over_nConstituents.size()/2];
328  }
329 
330  // debug info
331  if (Verbosity() > 0)
332  {
333  std::cout << "JetMultSub::EstimateRho(PHCompositeNode *topNode) Median et over nConstituents: " << median_et_over_nConstituents << std::endl;
334  }
335 
336  // return median et over nConstituents
337  return median_et_over_nConstituents;
338 }
339 
340 //____________________________________________________________________________..
341 // get multiplicity correction. Gets the multiplicity correction for a given pt
343 {
344  if (Verbosity() > 0)
345  {
346  std::cout << "JetMultSub::GetMultiplicityCorrection(float pt) Getting multiplicity correction" << std::endl;
347  }
348 
349 
350  // to do! Match jets and get the values for this function!
351  std::vector<float> pt_reco_bins = {0.0,100.0};
352  std::vector<float> multiplicity_corrections = {0.0};
353  // get pt bin
354  int pt_bin = -1;
355  int pt_size = pt_reco_bins.size();
356  for (int i = 0; i < pt_size; i++){
357  if (pt > pt_reco_bins[i]){
358  pt_bin = i;
359  }
360  }
361  float correction;
362  if(pt_bin == -1){
363  correction = 0.0;
364  }
365  else{
366  correction = multiplicity_corrections[pt_bin];
367  }
368  return correction;
369 }