Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationUpsilon.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationUpsilon.cc
1 #include "QAG4SimulationUpsilon.h"
2 
3 #include <qautils/QAHistManagerDef.h>
4 
5 #include <g4eval/SvtxEvalStack.h>
6 
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4TruthInfoContainer.h>
9 
10 #include <trackbase_historic/SvtxTrack.h>
11 
12 #include <g4eval/SvtxTrackEval.h> // for SvtxTrackEval
13 #include <g4eval/SvtxTruthEval.h> // for SvtxTruthEval
14 
15 #include <fun4all/Fun4AllHistoManager.h>
16 #include <fun4all/Fun4AllReturnCodes.h>
17 #include <fun4all/SubsysReco.h>
18 
19 #include <phool/getClass.h>
20 
21 #include <TAxis.h>
22 #include <TDatabasePDG.h>
23 #include <TH1.h>
24 #include <TH2.h>
25 #include <TNamed.h>
26 #include <TParticlePDG.h> // for TParticlePDG
27 #include <TString.h>
28 #include <TVector3.h>
29 
30 #include <CLHEP/Vector/LorentzVector.h>
31 #include <CLHEP/Vector/ThreeVector.h> // for Hep3Vector
32 
33 #include <cassert>
34 #include <cmath>
35 #include <cstdlib> // for abs
36 #include <iostream>
37 #include <utility> // for pair
38 
39 QAG4SimulationUpsilon::QAG4SimulationUpsilon(const std::string &name)
40  : SubsysReco(name)
41  , _svtxEvalStack(nullptr)
42  , m_etaRange(-1, 1)
43  , _truthContainer(nullptr)
44 {
45 }
46 
47 int QAG4SimulationUpsilon::InitRun(PHCompositeNode *topNode)
48 {
49  _truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
50  "G4TruthInfo");
51  if (!_truthContainer)
52  {
53  std::cout << "QAG4SimulationUpsilon::InitRun - Fatal Error - "
54  << "unable to find DST node "
55  << "G4TruthInfo" << std::endl;
56  assert(_truthContainer);
57  }
58 
59  if (!_svtxEvalStack)
60  {
61  _svtxEvalStack.reset(new SvtxEvalStack(topNode));
62  _svtxEvalStack->set_strict(true);
63  _svtxEvalStack->set_verbosity(Verbosity() + 1);
64  }
65 
66  return Fun4AllReturnCodes::EVENT_OK;
67 }
68 
69 int QAG4SimulationUpsilon::Init(PHCompositeNode * /*topNode*/)
70 {
71  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
72  assert(hm);
73 
74  // reco pT / gen pT histogram
75  TH1 *h(nullptr);
76 
77  h = new TH1F(TString(get_histo_prefix()) + "pTRecoGenRatio",
78  ";Reco p_{T}/Truth p_{T}", 500, 0, 2);
79  hm->registerHisto(h);
80 
81  h = new TH2F(TString(get_histo_prefix()) + "pTRecoGenRatio_pTGen",
82  ";Truth p_{T} [GeV/c];Reco p_{T}/Truth p_{T}", 200, 0, 50, 500, 0, 2);
83  // QAHistManagerDef::useLogBins(h->GetXaxis());
84  hm->registerHisto(h);
85 
86  // reco pT histogram
87  h = new TH1F(TString(get_histo_prefix()) + "nReco_pTGen",
88  "Reco tracks at truth p_{T};Truth p_{T} [GeV/c]", 200, 0.1, 50.5);
89  QAHistManagerDef::useLogBins(h->GetXaxis());
90  hm->registerHisto(h);
91  // reco pT histogram
92  h = new TH1F(TString(get_histo_prefix()) + "nGen_pTGen",
93  ";Truth p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
94  QAHistManagerDef::useLogBins(h->GetXaxis());
95  hm->registerHisto(h);
96 
97  // reco pT histogram
98  h = new TH1F(TString(get_histo_prefix()) + "nReco_etaGen",
99  "Reco tracks at truth #eta;Truth #eta", 200, -2, 2);
100  // QAHistManagerDef::useLogBins(h->GetXaxis());
101  hm->registerHisto(h);
102  // reco pT histogram
103  h = new TH1F(TString(get_histo_prefix()) + "nGen_etaGen",
104  ";Truth #eta;Track count / bin", 200, -2, 2);
105  // QAHistManagerDef::useLogBins(h->GetXaxis());
106  hm->registerHisto(h);
107 
108  h = new TH1F(TString(get_histo_prefix()) + "nGen_Pair_InvMassGen",
109  ";Truth Invariant Mass [GeV/c^2];Pair count / bin", 450, 0, 15);
110  // QAHistManagerDef::useLogBins(h->GetXaxis());
111  hm->registerHisto(h);
112  h = new TH1F(TString(get_histo_prefix()) + "nReco_Pair_InvMassReco",
113  ";Reco Invariant Mass [GeV/c^2];Pair count / bin", 450, 0, 15);
114  // QAHistManagerDef::useLogBins(h->GetXaxis());
115  hm->registerHisto(h);
116 
117  // n events and n tracks histogram
118  h = new TH1F(TString(get_histo_prefix()) + "Normalization",
119  TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
120  int i = 1;
121  h->GetXaxis()->SetBinLabel(i++, "Event");
122  h->GetXaxis()->SetBinLabel(i++, "Truth Track+");
123  h->GetXaxis()->SetBinLabel(i++, "Truth Track-");
124  h->GetXaxis()->SetBinLabel(i++, "Reco Track");
125  h->GetXaxis()->SetBinLabel(i++, "Truth Upsilon");
126  h->GetXaxis()->SetBinLabel(i++, "Truth Upsilon in Acc.");
127  h->GetXaxis()->SetBinLabel(i++, "Reco Upsilon");
128  h->GetXaxis()->LabelsOption("v");
129  hm->registerHisto(h);
130 
131  return Fun4AllReturnCodes::EVENT_OK;
132 }
133 
134 void QAG4SimulationUpsilon::addEmbeddingID(int embeddingID)
135 {
136  m_embeddingIDs.insert(embeddingID);
137 }
138 
139 int QAG4SimulationUpsilon::process_event(PHCompositeNode *topNode)
140 {
141  if (Verbosity() > 2)
142  std::cout << "QAG4SimulationUpsilon::process_event() entered" << std::endl;
143 
144  // histogram manager
145  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
146  assert(hm);
147 
148  if (_svtxEvalStack)
149  _svtxEvalStack->next_event(topNode);
150 
151  SvtxTrackEval *trackeval = _svtxEvalStack->get_track_eval();
152  assert(trackeval);
153  SvtxTruthEval *trutheval = _svtxEvalStack->get_truth_eval();
154  assert(trutheval);
155 
156  // reco pT / gen pT histogram
157  TH1 *h_pTRecoGenRatio = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio"));
158  assert(h_pTRecoGenRatio);
159 
160  // reco pT / gen pT histogram
161  TH2 *h_pTRecoGenRatio_pTGen = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio_pTGen"));
162  assert(h_pTRecoGenRatio);
163 
164  // reco histogram plotted at gen pT
165  TH1 *h_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTGen"));
166  assert(h_nReco_pTGen);
167 
168  // gen pT histogram
169  TH1 *h_nGen_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTGen"));
170  assert(h_nGen_pTGen);
171 
172  // reco histogram plotted at gen eta
173  TH1 *h_nReco_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_etaGen"));
174  assert(h_nReco_etaGen);
175 
176  // gen eta histogram
177  TH1 *h_nGen_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_etaGen"));
178  assert(h_nGen_etaGen);
179 
180  // inv mass
181  TH1 *h_nGen_Pair_InvMassGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_Pair_InvMassGen"));
182  assert(h_nGen_etaGen);
183  // inv mass
184  TH1 *h_nReco_Pair_InvMassReco = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_Pair_InvMassReco"));
185  assert(h_nGen_etaGen);
186 
187  // n events and n tracks histogram
188  TH1 *h_norm = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "Normalization"));
189  assert(h_norm);
190  h_norm->Fill("Event", 1);
191 
192  // buffer for daugther particles
193  typedef std::set<std::pair<PHG4Particle *, SvtxTrack *>> truth_reco_set_t;
194  truth_reco_set_t truth_reco_set_pos;
195  truth_reco_set_t truth_reco_set_neg;
196 
197  // fill histograms that need truth information
198  if (!_truthContainer)
199  {
200  std::cout << "QAG4SimulationUpsilon::process_event - fatal error - missing _truthContainer! ";
201  return Fun4AllReturnCodes::ABORTRUN;
202  }
203 
204  PHG4TruthInfoContainer::ConstRange range = _truthContainer->GetPrimaryParticleRange();
205  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
206  {
207  // get the truth particle information
208  PHG4Particle *g4particle = iter->second;
209 
210  if (Verbosity())
211  {
212  std::cout << "QAG4SimulationUpsilon::process_event - processing ";
213  g4particle->identify();
214  }
215 
216  if (m_embeddingIDs.size() > 0)
217  {
218  // only analyze subset of particle with proper embedding IDs
219  int candidate_embedding_id = trutheval->get_embed(g4particle);
220  if (candidate_embedding_id < 0) candidate_embedding_id = -1;
221 
222  // skip if no match
223  if (m_embeddingIDs.find(candidate_embedding_id) == m_embeddingIDs.end()) continue;
224  }
225 
226  double gpx = g4particle->get_px();
227  double gpy = g4particle->get_py();
228  double gpz = g4particle->get_pz();
229  double gpt = 0;
230  double geta = NAN;
231 
232  if (gpx != 0 && gpy != 0)
233  {
234  TVector3 gv(gpx, gpy, gpz);
235  gpt = gv.Pt();
236  geta = gv.Eta();
237  // gphi = gv.Phi();
238  }
239  if (m_etaRange.first < geta and geta < m_etaRange.second)
240  {
241  if (Verbosity())
242  {
243  std::cout << "QAG4SimulationUpsilon::process_event - accept particle eta = " << geta << std::endl;
244  }
245  }
246  else
247  {
248  if (Verbosity())
249  std::cout << "QAG4SimulationUpsilon::process_event - ignore particle eta = " << geta << std::endl;
250  continue;
251  }
252 
253  const int pid = g4particle->get_pid();
254 
255  if (abs(pid) != abs(m_daughterAbsPID))
256  {
257  if (Verbosity())
258  std::cout << "QAG4SimulationUpsilon::process_event - ignore particle PID = " << pid << " as m_daughterAbsPID = " << m_daughterAbsPID << std::endl;
259  continue;
260  }
261 
262  TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(pid);
263  if (!pdg_p)
264  {
265  std::cout << "QAG4SimulationUpsilon::process_event - Error - invalid particle ID = " << pid << std::endl;
266  continue;
267  }
268 
269  const double gcharge = pdg_p->Charge() / 3;
270  if (gcharge > 0)
271  {
272  h_norm->Fill("Truth Track+", 1);
273  }
274  else if (gcharge < 0)
275  {
276  h_norm->Fill("Truth Track-", 1);
277  }
278  else
279  {
280  if (Verbosity())
281  std::cout << "QAG4SimulationUpsilon::process_event - invalid neutral decay particle ID = " << pid << std::endl;
282  continue;
283  }
284  // h_norm->Fill("Truth Track", 1);
285 
286  h_nGen_pTGen->Fill(gpt);
287  h_nGen_etaGen->Fill(geta);
288 
289  // look for best matching track in reco data & get its information
290  SvtxTrack *track = trackeval->best_track_from(g4particle);
291  if (track)
292  {
293  h_nReco_etaGen->Fill(geta);
294  h_nReco_pTGen->Fill(gpt);
295 
296  double px = track->get_px();
297  double py = track->get_py();
298  double pz = track->get_pz();
299  double pt;
300  TVector3 v(px, py, pz);
301  pt = v.Pt();
302  // eta = v.Pt();
303  // phi = v.Pt();
304 
305  float pt_ratio = (gpt != 0) ? pt / gpt : 0;
306  h_pTRecoGenRatio->Fill(pt_ratio);
307  h_pTRecoGenRatio_pTGen->Fill(gpt, pt_ratio);
308  h_norm->Fill("Reco Track", 1);
309  }
310 
311  if (gcharge > 0)
312  {
313  truth_reco_set_pos.insert(std::make_pair(g4particle, track));
314  }
315  else if (gcharge < 0)
316  {
317  truth_reco_set_neg.insert(std::make_pair(g4particle, track));
318  }
319  } // for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
320 
321  // building pairs with buffer for daugter particles
322  TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(m_quarkoniaPID);
323  if (!pdg_p)
324  {
325  std::cout << "QAG4SimulationUpsilon::process_event - Fatal Error - invalid particle ID m_quarkoniaPID = " << m_quarkoniaPID << std::endl;
326 
327  return Fun4AllReturnCodes::ABORTRUN;
328  }
329  const double quarkonium_mass = pdg_p->Mass();
330  TParticlePDG *pdg_d = TDatabasePDG::Instance()->GetParticle(m_daughterAbsPID);
331  if (!pdg_d)
332  {
333  std::cout << "QAG4SimulationUpsilon::process_event - Fatal Error - invalid particle ID m_daughterAbsPID = " << m_daughterAbsPID << std::endl;
334 
335  return Fun4AllReturnCodes::ABORTRUN;
336  }
337  const double daughter_mass = pdg_d->Mass();
338 
339  for (const auto &pair_pos : truth_reco_set_pos)
340  for (const auto &pair_neg : truth_reco_set_neg)
341  {
342  assert(pair_pos.first);
343  assert(pair_neg.first);
344 
345  const CLHEP::HepLorentzVector gv_pos(
346  pair_pos.first->get_px(),
347  pair_pos.first->get_py(),
348  pair_pos.first->get_pz(),
349  pair_pos.first->get_e());
350 
351  const CLHEP::HepLorentzVector gv_neg(
352  pair_neg.first->get_px(),
353  pair_neg.first->get_py(),
354  pair_neg.first->get_pz(),
355  pair_neg.first->get_e());
356 
357  const CLHEP::HepLorentzVector gv_quakonium = gv_pos + gv_neg;
358 
359  if (fabs(quarkonium_mass - gv_quakonium.m()) > 1e-3)
360  {
361  if (Verbosity())
362  {
363  std::cout << "QAG4SimulationUpsilon::process_event - invalid pair with in compativle mass with " << quarkonium_mass << "GeV for PID = " << m_quarkoniaPID << ": " << std::endl;
364  pair_pos.first->identify();
365  pair_neg.first->identify();
366  }
367  continue;
368  }
369 
370  h_nGen_Pair_InvMassGen->Fill(gv_quakonium.m());
371  h_norm->Fill("Truth Upsilon in Acc.", 1);
372 
373  if (pair_pos.second and pair_neg.second)
374  {
375  CLHEP::HepLorentzVector v_pos;
376  CLHEP::HepLorentzVector v_neg;
377 
378  v_pos.setVectM(
379  CLHEP::Hep3Vector(
380  pair_pos.second->get_px(),
381  pair_pos.second->get_py(),
382  pair_pos.second->get_pz()),
383  daughter_mass);
384  v_neg.setVectM(
385  CLHEP::Hep3Vector(
386  pair_neg.second->get_px(),
387  pair_neg.second->get_py(),
388  pair_neg.second->get_pz()),
389  daughter_mass);
390 
391  const CLHEP::HepLorentzVector v_quakonium = v_pos + v_neg;
392 
393  h_nReco_Pair_InvMassReco->Fill(v_quakonium.m());
394  h_norm->Fill("Reco Upsilon", 1);
395 
396  } // if (pair_pos.second and pair_neg.second)
397 
398  } // for (const auto &pair_neg : truth_reco_set_neg)
399 
400  return Fun4AllReturnCodes::EVENT_OK;
401 }
402 
403 std::string
404 QAG4SimulationUpsilon::get_histo_prefix()
405 {
406  return std::string("h_") + Name() + std::string("_");
407 }