ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JetHepMCLoader.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetHepMCLoader.cc
1 // $Id: $
2 
11 #include "JetHepMCLoader.h"
12 
13 #include "JetMap.h" // for JetMap
14 #include "JetMapv1.h"
15 #include "Jetv1.h"
16 
19 
20 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY_A_LOT
23 #include <fun4all/Fun4AllServer.h>
24 #include <fun4all/SubsysReco.h> // for SubsysReco
25 
26 #include <phool/PHCompositeNode.h>
27 #include <phool/PHIODataNode.h>
28 #include <phool/PHNode.h> // for PHNode
29 #include <phool/PHNodeIterator.h>
30 #include <phool/PHObject.h> // for PHObject
31 #include <phool/getClass.h>
32 #include <phool/phool.h> // for PHWHERE
33 
34 #include <TAxis.h> // for TAxis
35 #include <TH1.h>
36 #include <TH2.h>
37 #include <TNamed.h> // for TNamed
38 
39 #include <HepMC/GenEvent.h> // for GenEvent, GenEvent::particl...
40 #include <HepMC/GenParticle.h> // for GenParticle
41 #include <HepMC/SimpleVector.h> // for FourVector
42 #include <HepMC/Units.h> // for conversion_factor, GEV
43 
44 #include <cassert>
45 #include <iostream>
46 
47 using namespace std;
48 
49 JetHepMCLoader::JetHepMCLoader(const std::string &jetInputCategory)
50  : SubsysReco("JetHepMCLoader_" + jetInputCategory)
51  , m_jetInputCategory(jetInputCategory)
52  , m_saveQAPlots(false)
53 {
54 }
55 
57 {
58 }
59 
61 {
63 }
64 
66 {
67  if (m_saveQAPlots)
68  {
70  assert(hm);
71 
72  const int n_bins = 1 + m_jetSrc.size();
73 
74  TH1D *h = new TH1D("hNormalization", //
75  "Normalization;Items;Summed quantity", n_bins, .5, n_bins + .5);
76  int i = 1;
77  h->GetXaxis()->SetBinLabel(i++, "Event count");
78  for (const hepmc_jet_src &src : m_jetSrc)
79  {
80  h->GetXaxis()->SetBinLabel(i++, (string("SubEvent ") + src.m_name).c_str());
81  }
82  h->GetXaxis()->LabelsOption("v");
83  hm->registerHisto(h);
84 
85  for (const hepmc_jet_src &src : m_jetSrc)
86  {
87  hm->registerHisto(
88  new TH2F((string("hJetEtEta_") + src.m_name).c_str(), //
89  ("Jet distribution from " + src.m_name + ";Jet #eta;Jet E_{T} [GeV]").c_str(),
90  40, -4., 4.,
91  100, 0, 100));
92  }
93  } // if (m_saveQAPlots)
94 
95  return CreateNodes(topNode);
96 }
97 
99 {
100  // For pile-up simulation: define GenEventMap
101  PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
102 
103  if (!genevtmap)
104  {
105  static bool once = true;
106 
107  if (once and Verbosity())
108  {
109  once = false;
110 
111  cout << "HepMCNodeReader::process_event - No PHHepMCGenEventMap node. Do not perform HepMC->Geant4 input" << endl;
112  }
113 
115  }
116 
117  if (m_saveQAPlots)
118  {
120  assert(hm);
121  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto("hNormalization"));
122  assert(h_norm);
123  h_norm->Fill("Event count", 1);
124  } // if (m_saveQAPlots)
125 
126  for (const hepmc_jet_src &src : m_jetSrc)
127  {
128  JetMap *jets = findNode::getClass<JetMap>(topNode, src.m_name);
129  assert(jets);
130 
131  jets->set_algo(src.m_algorithmID);
132  jets->set_par(src.m_parameter);
134 
135  PHHepMCGenEvent *genevt =
136  genevtmap->get(src.m_embeddingID);
137 
138  if (genevt == nullptr) continue;
139 
140  HepMC::GenEvent *evt = genevt->getEvent();
141  if (!evt)
142  {
143  cout << PHWHERE << " no evt pointer under HEPMC Node found:";
144  genevt->identify();
146  }
147 
148  assert(genevt->get_embedding_id() == src.m_embeddingID);
149 
150  TH2F *hjet = nullptr;
151 
152  if (m_saveQAPlots)
153  {
155  assert(hm);
156  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto("hNormalization"));
157  assert(h_norm);
158  h_norm->Fill((string("SubEvent ") + src.m_name).c_str(), 1);
159 
160  hjet = dynamic_cast<TH2F *>(hm->getHisto(string("hJetEtEta_") + src.m_name));
161  assert(hjet);
162 
163  } // if (m_saveQAPlots)
164 
165  const double mom_factor = HepMC::Units::conversion_factor(evt->momentum_unit(), HepMC::Units::GEV);
166 
167  for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
168  p != evt->particles_end(); ++p)
169  {
170  HepMC::GenParticle *part = (*p);
171 
172  assert(part);
173  if (Verbosity() >= VERBOSITY_A_LOT)
174  part->print();
175 
176  if (part->status() == src.m_tagStatus and part->pdg_id() == src.m_tagPID)
177  {
178  Jet *jet = new Jetv1();
179 
180  jet->set_px(part->momentum().px() * mom_factor);
181  jet->set_py(part->momentum().py() * mom_factor);
182  jet->set_pz(part->momentum().pz() * mom_factor);
183  jet->set_e(part->momentum().e() * mom_factor);
184 
185  jet->insert_comp(Jet::HEPMC_IMPORT, part->barcode());
186 
187  jets->insert(jet);
188 
189  if (hjet)
190  {
191  hjet->Fill(jet->get_eta(), jet->get_et());
192  }
193 
194  } // if (part->status() == src.m_tagStatus and part->pdg_id() == src.m_tagPID)
195  }
196  } // for (const hepmc_jet_src &src : m_jetSrc)
197 
199 }
200 
202 {
203  if (m_saveQAPlots)
204  {
206  assert(hm);
207 
208  cout << "JetHepMCLoader::End - saving QA histograms to " << Name() + ".root" << endl;
209  hm->dumpHistos(Name() + ".root", "RECREATE");
210  }
211 
213 }
214 
216  const std::string &name,
217  int embeddingID,
219  double parameter,
220  int tagPID,
221  int tagStatus)
222 {
223  string algorithmName = "Undefined_Jet_Algorithm";
224 
225  switch (algorithm)
226  {
227  case Jet::ANTIKT:
228  algorithmName = "ANTIKT";
229  break;
230 
231  case Jet::KT:
232  algorithmName = "KT";
233  break;
234 
235  case Jet::CAMBRIDGE:
236  algorithmName = "CAMBRIDGE";
237  break;
238 
239  default:
240 
241  break;
242  }
243 
244  hepmc_jet_src src{name, embeddingID, algorithmName, algorithm, parameter, tagPID, tagStatus};
245 
246  m_jetSrc.push_back(src);
247 }
248 
250 {
251  PHNodeIterator iter(topNode);
252 
253  // Looking for the DST node
254  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
255  if (!dstNode)
256  {
257  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
259  }
260 
261  for (const hepmc_jet_src &src : m_jetSrc)
262  {
263  // Create the AntiKt node if required
264  PHCompositeNode *AlgoNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", src.m_algorithmName.c_str()));
265  if (!AlgoNode)
266  {
267  AlgoNode = new PHCompositeNode(src.m_algorithmName.c_str());
268  dstNode->addNode(AlgoNode);
269  }
270 
271  // Create the Input node if required
272  PHCompositeNode *InputNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", m_jetInputCategory.c_str()));
273  if (!InputNode)
274  {
275  InputNode = new PHCompositeNode(m_jetInputCategory.c_str());
276  AlgoNode->addNode(InputNode);
277  }
278 
279  JetMap *jets = findNode::getClass<JetMap>(topNode, src.m_name);
280  if (!jets)
281  {
282  jets = new JetMapv1();
283  PHIODataNode<PHObject> *JetMapNode = new PHIODataNode<PHObject>(jets, src.m_name.c_str(), "PHObject");
284  InputNode->addNode(JetMapNode);
285  }
286  }
287 
289 }
290 
293 {
294  string histname(Name() + "_HISTOS");
295 
297  Fun4AllHistoManager *hm = se->getHistoManager(histname);
298 
299  if (not hm)
300  {
301  cout
302  << "TPCDataStreamEmulator::get_HistoManager - Making Fun4AllHistoManager " << histname
303  << endl;
304  hm = new Fun4AllHistoManager(histname);
305  se->registerHistoManager(hm);
306  }
307 
308  assert(hm);
309 
310  return hm;
311 }