ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFParticle_sPHENIX.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFParticle_sPHENIX.cc
1 /*
2  * This file is part of KFParticle package
3  * Copyright ( C ) 2007-2019 FIAS Frankfurt Institute for Advanced Studies
4  * 2007-2019 Goethe University of Frankfurt
5  * 2007-2019 Ivan Kisel <I.Kisel@compeng.uni-frankfurt.de>
6  * 2007-2019 Maksym Zyzak
7  *
8  * KFParticle is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * ( at your option ) any later version.
12  *
13  * KFParticle is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program. If not, see <https://www.gnu.org/licenses/>.
20  */
21 
22 #include "KFParticle_sPHENIX.h"
23 
26 
28 
29 #include <phool/getClass.h>
30 
31 #include <TFile.h>
32 
33 typedef std::pair<int, float> particle_pair;
34 
37 std::map<std::string, particle_pair> particleList = kfp_list.getParticleList();
39 
42  : SubsysReco("KFPARTICLE")
43  , m_has_intermediates_sPHENIX(false)
44  , m_constrain_to_vertex_sPHENIX(false)
45  , m_require_mva(false)
46  , m_save_dst(0)
47  , m_save_output(1)
48  , m_outfile_name("outputData.root")
49  , m_outfile(nullptr)
50 {
51 }
52 
54  : SubsysReco(name)
55  , m_has_intermediates_sPHENIX(false)
56  , m_constrain_to_vertex_sPHENIX(false)
57  , m_require_mva(false)
58  , m_save_dst(0)
59  , m_save_output(1)
60  , m_outfile_name("outputData.root")
61  , m_outfile(nullptr)
62 {
63 }
64 
66 {
67  if (m_save_output && Verbosity() >= VERBOSITY_SOME) std::cout << "Output nTuple: " << m_outfile_name << std::endl;
68 
69  if (m_save_dst) createParticleNode(topNode);
70 
71  if (m_require_mva)
72  {
73  TMVA::Reader *reader;
74  std::vector<Float_t> MVA_parValues;
75  tie(reader, MVA_parValues) = initMVA();
76  }
77 
78  for (int i = 0; i < m_num_tracks; ++i)
79  if (!particleList.count(m_daughter_name[i]))
80  {
81  std::cout << "Your track PID, " << m_daughter_name[i] << "is not in the particle list" << std::endl;
82  std::cout << "Check KFParticle_particleList.cxx for a list of available particles" << std::endl;
83  exit(0);
84  }
85 
86  int returnCode = 0;
87  if (!m_decayDescriptor.empty()) returnCode = parseDecayDescriptor();
88 
89  return returnCode;
90 }
91 
93 {
94  std::vector<KFParticle> mother, vertex;
95  std::vector<std::vector<KFParticle>> daughters, intermediates;
96  int nPVs, multiplicity;
97 
98  if (!m_use_fake_pv)
99  {
100  SvtxVertexMap *check_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vtx_map_node_name);
101  if (check_vertexmap->size() == 0)
102  {
103  if (Verbosity() >= VERBOSITY_SOME) std::cout << "KFParticle: Event skipped as there are no vertices" << std::endl;
105  }
106  }
107 
108  SvtxTrackMap *check_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name);
109  if (check_trackmap->size() == 0)
110  {
111  if (Verbosity() >= VERBOSITY_SOME) std::cout << "KFParticle: Event skipped as there are no tracks" << std::endl;
113  }
114 
115  createDecay(topNode, mother, vertex, daughters, intermediates, nPVs, multiplicity);
116 
117  if (!m_has_intermediates_sPHENIX) intermediates = daughters;
118  if (!m_constrain_to_vertex_sPHENIX) vertex = mother;
119 
120  if (mother.size() != 0)
121  {
122  for (unsigned int i = 0; i < mother.size(); ++i)
123  {
124  if (m_save_output && candidateCounter == 0)
125  {
126  m_outfile = new TFile(m_outfile_name.c_str(), "RECREATE");
128  }
129 
130  candidateCounter += 1;
131 
132  if (m_save_output) fillBranch(topNode, mother[i], vertex[i], daughters[i], intermediates[i], nPVs, multiplicity);
133  if (m_save_dst) fillParticleNode(topNode, mother[i], daughters[i], intermediates[i]);
134 
135  if (Verbosity() >= VERBOSITY_SOME)
136  {
137  printParticles(mother[i], vertex[i], daughters[i], intermediates[i], nPVs, multiplicity);
138  }
139  if (Verbosity() >= VERBOSITY_MORE)
140  {
141  if (m_save_dst) printNode(topNode);
142  }
143  }
144  }
145 
147 }
148 
150 {
151  std::cout << "KFParticle_sPHENIX object " << Name() << " finished. Number of canadidates: " << candidateCounter << std::endl;
152 
153  if (m_save_output && candidateCounter != 0)
154  {
155  m_outfile->Write();
156  m_outfile->Close();
157  delete m_outfile;
158  }
159 
160  return 0;
161 }
162 
163 void KFParticle_sPHENIX::printParticles(KFParticle motherParticle,
164  KFParticle chosenVertex,
165  std::vector<KFParticle> daughterParticles,
166  std::vector<KFParticle> intermediateParticles,
167  int numPVs, int numTracks)
168 {
169  std::cout << "\n---------------KFParticle candidate information---------------" << std::endl;
170 
171  std::cout << "Mother information:" << std::endl;
172  identify(motherParticle);
173 
175  {
176  std::cout << "Intermediate state information:" << std::endl;
177  for (unsigned int i = 0; i < intermediateParticles.size(); i++)
178  {
179  identify(intermediateParticles[i]);
180  }
181  }
182 
183  std::cout << "Final track information:" << std::endl;
184  for (unsigned int i = 0; i < daughterParticles.size(); i++)
185  {
186  identify(daughterParticles[i]);
187  }
188 
190  {
191  std::cout << "Primary vertex information:" << std::endl;
192  std::cout << "(x,y,z) = (" << chosenVertex.GetX() << " +/- " << sqrt(chosenVertex.GetCovariance(0, 0)) << ", ";
193  std::cout << chosenVertex.GetY() << " +/- " << sqrt(chosenVertex.GetCovariance(1, 1)) << ", ";
194  std::cout << chosenVertex.GetZ() << " +/- " << sqrt(chosenVertex.GetCovariance(2, 2)) << ") cm\n"
195  << std::endl;
196  }
197 
198  std::cout << "The number of primary vertices is: " << numPVs << std::endl;
199  std::cout << "The number of tracks in the event is: " << numTracks << std::endl;
200 
201  std::cout << "------------------------------------------------------------\n"
202  << std::endl;
203 }
204 
206 {
207  bool ddCanBeParsed = true;
208 
209  size_t daughterLocator;
210 
211  std::string mother;
212  std::string intermediate;
213  std::string daughter;
214 
215  std::vector<std::pair<std::string, int>> intermediate_list;
216  std::vector<std::string> intermediates_name;
217  std::vector<int> intermediates_charge;
218 
219  std::vector<std::pair<std::string, int>> daughter_list;
220  std::vector<std::string> daughters_name;
221  std::vector<int> daughters_charge;
222 
223  int nTracks = 0;
224  std::vector<int> m_nTracksFromIntermediates;
225 
226  std::string decayArrow = "->";
227  std::string chargeIndicator = "^";
228  std::string startIntermediate = "{";
229  std::string endIntermediate = "}";
230 
231  std::string manipulateDecayDescriptor = m_decayDescriptor;
232 
233  //Remove all white space before we begin
234  size_t pos;
235  while ((pos = manipulateDecayDescriptor.find(" ")) != std::string::npos) manipulateDecayDescriptor.replace(pos, 1, "");
236 
237  //Check for charge conjugate requirement
238  std::string checkForCC = manipulateDecayDescriptor.substr(0, 1) + manipulateDecayDescriptor.substr(manipulateDecayDescriptor.size()-3, 3);
239  std::for_each(checkForCC.begin(), checkForCC.end(), [](char & c) {c = ::toupper(c);});
240 
241  //Remove the CC check if needed
242  if (checkForCC == "[]CC")
243  {
244  manipulateDecayDescriptor = manipulateDecayDescriptor.substr(1, manipulateDecayDescriptor.size()-4);
245  getChargeConjugate(true);
246  }
247 
248  //Find the initial particle
249  size_t findMotherEndPoint = manipulateDecayDescriptor.find(decayArrow);
250  mother = manipulateDecayDescriptor.substr(0, findMotherEndPoint);
251  if (!findParticle(mother)) ddCanBeParsed = false;
252  manipulateDecayDescriptor.erase(0, findMotherEndPoint + decayArrow.length());
253 
254  //Try and find the intermediates
255  while ((pos = manipulateDecayDescriptor.find(startIntermediate)) != std::string::npos)
256  {
257  size_t findIntermediateStartPoint = manipulateDecayDescriptor.find(startIntermediate, pos);
258  size_t findIntermediateEndPoint = manipulateDecayDescriptor.find(endIntermediate, pos);
259  std::string intermediateDecay = manipulateDecayDescriptor.substr(pos + 1, findIntermediateEndPoint - (pos + 1));
260 
261  intermediate = intermediateDecay.substr(0, intermediateDecay.find(decayArrow));
262  if (findParticle(intermediate)) intermediates_name.push_back(intermediate.c_str());
263  else ddCanBeParsed = false;
264 
265  //Now find the daughters associated to this intermediate
266  int nDaughters = 0;
267  intermediateDecay.erase(0, intermediateDecay.find(decayArrow) + decayArrow.length());
268  while ((daughterLocator = intermediateDecay.find(chargeIndicator)) != std::string::npos)
269  {
270  daughter = intermediateDecay.substr(0, daughterLocator);
271  if (findParticle(daughter))
272  {
273  daughters_name.push_back(daughter.c_str());
274 
275  std::string daughterChargeString = intermediateDecay.substr(daughterLocator + 1, 1);
276  if (daughterChargeString == "+")
277  {
278  daughters_charge.push_back(+1);
279  }
280  else if (daughterChargeString == "-")
281  {
282  daughters_charge.push_back(-1);
283  }
284  else if (daughterChargeString == "0")
285  {
286  daughters_charge.push_back(0);
287  }
288  else
289  {
290  if (Verbosity() >= VERBOSITY_MORE) std::cout << "The charge of " << daughterChargeString << " was not known" << std::endl;
291  ddCanBeParsed = false;
292  }
293  }
294  else ddCanBeParsed = false;
295  intermediateDecay.erase(0, daughterLocator + 2);
296  ++nDaughters;
297  }
298  manipulateDecayDescriptor.erase(findIntermediateStartPoint, findIntermediateEndPoint + 1 - findIntermediateStartPoint);
299  m_nTracksFromIntermediates.push_back(nDaughters);
300  nTracks += nDaughters;
301  }
302 
303  //Now find any remaining reconstructable tracks from the mother
304  while ((daughterLocator = manipulateDecayDescriptor.find(chargeIndicator)) != std::string::npos)
305  {
306  daughter = manipulateDecayDescriptor.substr(0, daughterLocator);
307  if (findParticle(daughter))
308  {
309  daughters_name.push_back(daughter.c_str());
310  std::string daughterChargeString = manipulateDecayDescriptor.substr(daughterLocator + 1, 1);
311  if (daughterChargeString == "+")
312  {
313  daughters_charge.push_back(+1);
314  }
315  else if (daughterChargeString == "-")
316  {
317  daughters_charge.push_back(-1);
318  }
319  else if (daughterChargeString == "0")
320  {
321  daughters_charge.push_back(0);
322  }
323  else
324  {
325  if (Verbosity() >= VERBOSITY_MORE) std::cout << "The charge of " << daughterChargeString << " was not known" << std::endl;
326  ddCanBeParsed = false;
327  }
328  }
329  else ddCanBeParsed = false;
330  manipulateDecayDescriptor.erase(0, daughterLocator + 2);
331  nTracks += 1;
332  }
333 
334  int trackStart = 0;
335  int trackEnd = 0;
336  for (unsigned int i = 0; i < intermediates_name.size(); ++i)
337  {
338  trackStart = trackEnd;
339  trackEnd = m_nTracksFromIntermediates[i] + trackStart;
340 
341  int vtxCharge = 0;
342 
343  for (int j = trackStart; j < trackEnd; ++j)
344  {
345  vtxCharge += daughters_charge[j];
346  }
347 
348  intermediates_charge.push_back(vtxCharge);
349 
350  intermediate_list.push_back(std::make_pair(intermediates_name[i], intermediates_charge[i]));
351  }
352 
353  for (int i = 0; i < nTracks; ++i)
354  {
355  daughter_list.push_back(std::make_pair(daughters_name[i], daughters_charge[i]));
356  }
357 
358  setMotherName(mother);
359  setNumberOfTracks(nTracks);
360  setDaughters(daughter_list);
361 
362  if (intermediates_name.size() > 0)
363  {
364  hasIntermediateStates(true);
365  setIntermediateStates(intermediate_list);
366  setNumberOfIntermediateStates(intermediates_name.size());
367  setNumberTracksFromIntermeditateState(m_nTracksFromIntermediates);
368  }
369 
370  if (ddCanBeParsed)
371  {
372  if (Verbosity() >= VERBOSITY_MORE) std::cout << "Your decay descriptor can be parsed" << std::endl;
373  return 0;
374  }
375  else
376  {
377  if (Verbosity() >= VERBOSITY_SOME) std::cout << "KFParticle: Your decay descriptor, " << Name() << " cannot be parsed" << "\nExiting!" << std::endl;
379  }
380 }
381 
383 {
384  bool particleFound = true;
385  if (!particleList.count(particle))
386  {
387  if (Verbosity() >= VERBOSITY_SOME)
388  {
389  std::cout << "The particle, " << particle << " is not recognized" << std::endl;
390  std::cout << "Check KFParticle_particleList.cc for a list of available particles" << std::endl;
391  }
392  particleFound = false;
393  }
394 
395  return particleFound;
396 }