ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JetReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetReco.cc
1 
2 #include "JetReco.h"
3 
4 #include "Jet.h"
5 #include "JetAlgo.h"
6 #include "JetInput.h"
7 #include "JetMap.h"
8 #include "JetMapv1.h"
9 
10 // PHENIX includes
12 #include <fun4all/SubsysReco.h> // for SubsysReco
13 
14 #include <phool/PHCompositeNode.h>
15 #include <phool/PHIODataNode.h>
16 #include <phool/PHNode.h> // for PHNode
17 #include <phool/PHNodeIterator.h>
19 #include <phool/PHObject.h> // for PHObject
20 #include <phool/getClass.h>
21 #include <phool/phool.h> // for PHWHERE
22 
23 // standard includes
24 #include <cstdlib> // for exit
25 #include <iostream>
26 #include <memory> // for allocator_traits<>::value_type
27 #include <vector>
28 
29 using namespace std;
30 
31 JetReco::JetReco(const string &name)
32  : SubsysReco(name)
33  , _inputs()
34  , _algos()
35  , _algonode()
36  , _inputnode()
37  , _outputs()
38 {
39 }
40 
42 {
43  for (unsigned int i = 0; i < _inputs.size(); ++i) delete _inputs[i];
44  _inputs.clear();
45  for (unsigned int i = 0; i < _algos.size(); ++i) delete _algos[i];
46  _algos.clear();
47  _outputs.clear();
48 }
49 
50 int JetReco::Init(PHCompositeNode */*topNode*/)
51 {
53 }
54 
56 {
57  if (Verbosity() > 0)
58  {
59  cout << "========================== JetReco::InitRun() =============================" << endl;
60  cout << " Input Selections:" << endl;
61  for (unsigned int i = 0; i < _inputs.size(); ++i) _inputs[i]->identify();
62  cout << " Algorithms:" << endl;
63  for (unsigned int i = 0; i < _algos.size(); ++i) _algos[i]->identify();
64  cout << "===========================================================================" << endl;
65  }
66 
67  return CreateNodes(topNode);
68 }
69 
71 {
72  if (Verbosity() > 1) cout << "JetReco::process_event -- entered" << endl;
73 
74  //---------------------------------
75  // Get Objects off of the Node Tree
76  //---------------------------------
77 
78  std::vector<Jet *> inputs; // owns memory
79  for (unsigned int iselect = 0; iselect < _inputs.size(); ++iselect)
80  {
81  std::vector<Jet *> parts = _inputs[iselect]->get_input(topNode);
82  for (unsigned int ipart = 0; ipart < parts.size(); ++ipart)
83  {
84  inputs.push_back(parts[ipart]);
85  inputs.back()->set_id(inputs.size() - 1); // unique ids ensured
86  }
87  }
88 
89  //---------------------------
90  // Run the jet reconstruction
91  //---------------------------
92  for (unsigned int ialgo = 0; ialgo < _algos.size(); ++ialgo)
93  {
94  std::vector<Jet *> jets = _algos[ialgo]->get_jets(inputs); // owns memory
95 
96  // send the output somewhere on the DST
97  FillJetNode(topNode, ialgo, jets);
98  }
99 
100  // clean up input vector
101  for (unsigned int i = 0; i < inputs.size(); ++i) delete inputs[i];
102  inputs.clear();
103 
104  if (Verbosity() > 1) cout << "JetReco::process_event -- exited" << endl;
105 
107 }
108 
109 int JetReco::End(PHCompositeNode */*topNode*/)
110 {
112 }
113 
115 {
116  PHNodeIterator iter(topNode);
117 
118  // Looking for the DST node
119  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
120  if (!dstNode)
121  {
122  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
124  }
125 
126  // Create the AntiKt node if required
127  PHCompositeNode *AlgoNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _algonode.c_str()));
128  if (!AlgoNode)
129  {
130  AlgoNode = new PHCompositeNode(_algonode.c_str());
131  dstNode->addNode(AlgoNode);
132  }
133 
134  // Create the Input node if required
135  PHCompositeNode *InputNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _inputnode.c_str()));
136  if (!InputNode)
137  {
138  InputNode = new PHCompositeNode(_inputnode.c_str());
139  AlgoNode->addNode(InputNode);
140  }
141 
142  for (unsigned i = 0; i < _outputs.size(); ++i)
143  {
144  JetMap *jets = findNode::getClass<JetMap>(topNode, _outputs[i]);
145  if (!jets)
146  {
147  jets = new JetMapv1();
148  PHIODataNode<PHObject> *JetMapNode = new PHIODataNode<PHObject>(jets, _outputs[i].c_str(), "PHObject");
149  InputNode->addNode(JetMapNode);
150  }
151  }
152 
154 }
155 
156 void JetReco::FillJetNode(PHCompositeNode *topNode, int ipos, std::vector<Jet *> jets)
157 {
158  JetMap *jetmap = nullptr;
159  PHTypedNodeIterator<JetMap> jetmapiter(topNode);
160  PHIODataNode<JetMap> *JetMapNode = jetmapiter.find(_outputs[ipos].c_str());
161  if (!JetMapNode)
162  {
163  cout << PHWHERE << " ERROR: Can't find JetMap: " << _outputs[ipos] << endl;
164  exit(-1);
165  }
166  else
167  {
168  jetmap = dynamic_cast<JetMap *> (JetMapNode->getData());
169  }
170 
171  jetmap->set_algo(_algos[ipos]->get_algo());
172  jetmap->set_par(_algos[ipos]->get_par());
173  for (unsigned int i = 0; i < _inputs.size(); ++i)
174  {
175  jetmap->insert_src(_inputs[i]->get_src());
176  }
177 
178  for (unsigned int i = 0; i < jets.size(); ++i)
179  {
180  jetmap->insert(jets[i]); // map takes ownership, sets unique id
181  }
182 
183  return;
184 }