ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KFParticle_Tools.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KFParticle_Tools.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 /*****************/
23 /* Cameron Dean */
24 /* LANL 2020 */
25 /* cdean@bnl.gov */
26 /*****************/
27 
28 #include "KFParticle_Tools.h"
29 
30 #include <phool/getClass.h>
35 
36 //KFParticle stuff
37 #include <KFPTrack.h>
38 #include <KFParticle.h>
39 #include <KFParticleDatabase.h>
40 #include <KFVertex.h>
41 
42 #include <TMatrixD.h>
43 
44 #include <Eigen/Core>
45 #include <Eigen/Dense>
46 
48 typedef std::pair<int, float> particle_pair;
50 
51 //Particle masses are in GeV
52 std::map<std::string, particle_pair> particleMasses = kfp_particleList.getParticleList();
53 
56  : m_has_intermediates(false)
57  , m_min_mass(0)
58  , m_max_mass(0)
59  , m_min_decayTime(-1*FLT_MAX)
60  , m_max_decayTime(FLT_MAX)
61  , m_min_decayLength(-1*FLT_MAX)
62  , m_max_decayLength(FLT_MAX)
63  , m_track_pt(0.25)
64  , m_track_ptchi2(FLT_MAX)
65  , m_track_ip(-1.)
66  , m_track_ipchi2(10.)
67  , m_track_chi2ndof(4.)
68  , m_comb_DCA(0.05)
69  , m_vertex_chi2ndof(15.)
70  , m_fdchi2(50.)
71  , m_dira_min(0.95)
72  , m_dira_max(1.01)
73  , m_mother_pt(0.)
74  , m_mother_ipchi2(FLT_MAX)
75  , m_get_charge_conjugate(true)
76  , m_vtx_map_node_name("SvtxVertexMap")
77  , m_trk_map_node_name("SvtxTrackMap")
78  , m_dst_vertexmap()
79  , m_dst_trackmap()
80  , m_dst_vertex()
81  , m_dst_track()
82 {
83 }
84 
86 {
87  float f_vertexParameters[6] = {m_dst_vertex->get_x(),
89  m_dst_vertex->get_z(), 0, 0, 0};
90 
91  float f_vertexCovariance[21];
92  unsigned int iterate = 0;
93  for (unsigned int i = 0; i < 3; ++i)
94  for (unsigned int j = 0; j <= i; ++j)
95  {
96  f_vertexCovariance[iterate] = m_dst_vertex->get_error(i, j);
97  ++iterate;
98  }
99 
100  KFParticle kfp_vertex;
101  kfp_vertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
102  kfp_vertex.NDF() = m_dst_vertex->get_ndof();
103  kfp_vertex.Chi2() = m_dst_vertex->get_chisq();
104 
105  return kfp_vertex;
106 }
107 
108 std::vector<KFParticle> KFParticle_Tools::makeAllPrimaryVertices(PHCompositeNode *topNode, std::string vertexMapName)
109 {
110  std::string vtxMN;
111  if (vertexMapName.empty())
112  vtxMN = m_vtx_map_node_name;
113  else
114  vtxMN = vertexMapName;
115 
116  std::vector<KFParticle> primaryVertices;
117  m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, vtxMN);
118  unsigned int vertexID = 0;
119 
120  for (SvtxVertexMap::ConstIter iter = m_dst_vertexmap->begin(); iter != m_dst_vertexmap->end(); ++iter)
121  {
122  m_dst_vertex = iter->second;
123 
124  primaryVertices.push_back(makeVertex(topNode));
125  primaryVertices[vertexID].SetId(iter->first);
126  ++vertexID;
127  }
128 
129  return primaryVertices;
130 }
131 
133 {
134  float f_trackParameters[6] = {m_dst_track->get_x(),
135  m_dst_track->get_y(),
136  m_dst_track->get_z(),
137  m_dst_track->get_px(),
138  m_dst_track->get_py(),
139  m_dst_track->get_pz()};
140 
141  float f_trackCovariance[21];
142  unsigned int iterate = 0;
143  for (unsigned int i = 0; i < 6; ++i)
144  for (unsigned int j = 0; j <= i; ++j)
145  {
146  f_trackCovariance[iterate] = m_dst_track->get_error(i, j);
147  ++iterate;
148  }
149 
150  KFParticle kfp_particle;
151  kfp_particle.Create(f_trackParameters, f_trackCovariance, (Int_t) m_dst_track->get_charge(), -1);
152  kfp_particle.NDF() = m_dst_track->get_ndf();
153  kfp_particle.Chi2() = m_dst_track->get_chisq();
154  kfp_particle.SetId(m_dst_track->get_id());
155 
156  return kfp_particle;
157 }
158 
160 {
161  std::vector<KFParticle> daughterParticles;
162  m_dst_trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trk_map_node_name.c_str());
163  unsigned int trackID = 0;
164 
165  for (SvtxTrackMap::Iter iter = m_dst_trackmap->begin(); iter != m_dst_trackmap->end(); ++iter)
166  {
167  m_dst_track = iter->second;
168  daughterParticles.push_back(makeParticle(topNode));
169  daughterParticles[trackID].SetId(iter->first);
170  ++trackID;
171  }
172 
173  return daughterParticles;
174 }
175 
176 int KFParticle_Tools::getTracksFromVertex(PHCompositeNode *topNode, KFParticle vertex, std::string vertexMapName)
177 {
178  std::string vtxMN;
179  if (vertexMapName.empty())
180  vtxMN = m_vtx_map_node_name;
181  else
182  vtxMN = vertexMapName;
183 
184  SvtxVertex *associatedVertex = NULL;
185  m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, vtxMN);
186 
187  associatedVertex = m_dst_vertexmap->find(vertex.Id())->second;
188 
189  return associatedVertex->size_tracks();
190 }
191 
192 /*const*/ bool KFParticle_Tools::isGoodTrack(KFParticle particle, const std::vector<KFParticle> primaryVertices)
193 {
194  bool goodTrack = false;
195 
196  float min_ip = 0;
197  float min_ipchi2 = 0;
198 
199  float pt = particle.GetPt();
200  float pterr = particle.GetErrPt();
201  float ptchi2 = pow(pterr / pt, 2);
202  float trackchi2ndof = particle.GetChi2() / particle.GetNDF();
203  calcMinIP(particle, primaryVertices, min_ip, min_ipchi2);
204 
205  if (pt >= m_track_pt && ptchi2 <= m_track_ptchi2
206  && min_ip >= m_track_ip && min_ipchi2 >= m_track_ipchi2
207  && trackchi2ndof <= m_track_chi2ndof)
208  goodTrack = true;
209 
210  return goodTrack;
211 }
212 
213 int KFParticle_Tools::calcMinIP(KFParticle track, std::vector<KFParticle> PVs,
214  float& minimumIP, float& minimumIPchi2)
215 {
216  std::vector<float> ip, ipchi2;
217 
218  for (unsigned int i_verts = 0; i_verts < PVs.size(); ++i_verts)
219  {
220  ip.push_back(track.GetDistanceFromVertex(PVs[i_verts]));
221  ipchi2.push_back(track.GetDeviationFromVertex(PVs[i_verts]));
222  }
223 
224  auto minmax_ip = minmax_element(ip.begin(), ip.end()); //Order the IP from small to large
225  minimumIP = *minmax_ip.first;
226  auto minmax_ipchi2 = minmax_element(ipchi2.begin(), ipchi2.end()); //Order the IP chi2 from small to large
227  minimumIPchi2 = *minmax_ipchi2.first;
228 
229  return 0;
230 }
231 
232 std::vector<int> KFParticle_Tools::findAllGoodTracks(std::vector<KFParticle> daughterParticles, const std::vector<KFParticle> primaryVertices)
233 {
234  std::vector<int> goodTrackIndex;
235 
236  for (unsigned int i_parts = 0; i_parts < daughterParticles.size(); ++i_parts)
237  {
238  if (isGoodTrack(daughterParticles[i_parts], primaryVertices)) goodTrackIndex.push_back(i_parts);
239  }
240 
241  removeDuplicates(goodTrackIndex);
242 
243  return goodTrackIndex;
244 }
245 
246 std::vector<std::vector<int>> KFParticle_Tools::findTwoProngs(std::vector<KFParticle> daughterParticles, std::vector<int> goodTrackIndex, int nTracks)
247 {
248  std::vector<std::vector<int>> goodTracksThatMeet;
249 
250  for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
251  {
252  for (std::vector<int>::iterator j_it = goodTrackIndex.begin(); j_it != goodTrackIndex.end(); ++j_it)
253  {
254  if (i_it < j_it)
255  {
256  if (daughterParticles[*i_it].GetDistanceFromParticle(daughterParticles[*j_it]) <= m_comb_DCA)
257  {
258  KFVertex twoParticleVertex;
259  twoParticleVertex += daughterParticles[*i_it];
260  twoParticleVertex += daughterParticles[*j_it];
261  float vertexchi2ndof = twoParticleVertex.GetChi2() / twoParticleVertex.GetNDF();
262  std::vector<int> combination = {*i_it, *j_it};
263 
264  if (nTracks == 2 && vertexchi2ndof <= m_vertex_chi2ndof)
265  {
266  goodTracksThatMeet.push_back(combination);
267  }
268  else if (nTracks == 2 && vertexchi2ndof > m_vertex_chi2ndof)
269  {
270  continue;
271  }
272  else
273  {
274  goodTracksThatMeet.push_back(combination);
275  }
276  }
277  }
278  }
279  }
280 
281  return goodTracksThatMeet;
282 }
283 
284 std::vector<std::vector<int>> KFParticle_Tools::findNProngs(std::vector<KFParticle> daughterParticles,
285  std::vector<int> goodTrackIndex,
286  std::vector<std::vector<int>> goodTracksThatMeet,
287  int nRequiredTracks, unsigned int nProngs)
288 {
289  unsigned int nGoodProngs = goodTracksThatMeet.size();
290 
291  for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
292  {
293  for (unsigned int i_prongs = 0; i_prongs < nGoodProngs; ++i_prongs)
294  {
295  bool trackNotUsedAlready = true;
296  for (unsigned int i_trackCheck = 0; i_trackCheck < nProngs - 1; ++i_trackCheck)
297  if (*i_it == goodTracksThatMeet[i_prongs][i_trackCheck]) trackNotUsedAlready = false;
298  if (trackNotUsedAlready)
299  {
300  bool dcaMet = 1;
301  for (unsigned int i = 0; i < nProngs - 1; ++i)
302  {
303  if (daughterParticles[*i_it].GetDistanceFromParticle(daughterParticles[goodTracksThatMeet[i_prongs][i]]) > m_comb_DCA)
304  {
305  dcaMet = 0;
306  }
307  }
308 
309  if (dcaMet)
310  {
311  KFVertex particleVertex;
312  particleVertex += daughterParticles[*i_it];
313  std::vector<int> combination;
314  combination.push_back(*i_it);
315  for (unsigned int i = 0; i < nProngs - 1; ++i)
316  {
317  particleVertex += daughterParticles[goodTracksThatMeet[i_prongs][i]];
318  combination.push_back(goodTracksThatMeet[i_prongs][i]);
319  }
320  float vertexchi2ndof = particleVertex.GetChi2() / particleVertex.GetNDF();
321 
322  if ((unsigned int) nRequiredTracks == nProngs && vertexchi2ndof <= m_vertex_chi2ndof)
323  {
324  goodTracksThatMeet.push_back(combination);
325  }
326  else if ((unsigned int) nRequiredTracks == nProngs && vertexchi2ndof > m_vertex_chi2ndof)
327  {
328  continue;
329  }
330  else
331  {
332  goodTracksThatMeet.push_back(combination);
333  }
334  }
335  }
336  }
337  }
338 
339  goodTracksThatMeet.erase(goodTracksThatMeet.begin(), goodTracksThatMeet.begin() + nGoodProngs);
340  for (unsigned int i = 0; i < goodTracksThatMeet.size(); ++i) sort(goodTracksThatMeet[i].begin(), goodTracksThatMeet[i].end());
341  removeDuplicates(goodTracksThatMeet);
342 
343  return goodTracksThatMeet;
344 }
345 
346 std::vector<std::vector<int>> KFParticle_Tools::appendTracksToIntermediates(KFParticle intermediateResonances[], std::vector<KFParticle> daughterParticles, std::vector<int> goodTrackIndex, int num_remaining_tracks)
347 {
348  std::vector<std::vector<int>> goodTracksThatMeet, goodTracksThatMeetIntermediates; //, vectorOfGoodTracks;
349  if (num_remaining_tracks == 1)
350  {
351  for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
352  {
353  std::vector<KFParticle> v_intermediateResonances(intermediateResonances, intermediateResonances + m_num_intermediate_states);
354  std::vector<std::vector<int>> dummyTrackList;
355  std::vector<int> dummyTrackID; //I already have the track ids stored in goodTracksThatMeet[i]
356  v_intermediateResonances.insert(end(v_intermediateResonances), daughterParticles[*i_it]);
357  for (unsigned int k = 0; k < v_intermediateResonances.size(); ++k)
358  {
359  dummyTrackID.push_back(k);
360  }
361  dummyTrackList = findTwoProngs(v_intermediateResonances, dummyTrackID, (int) v_intermediateResonances.size());
362  if (v_intermediateResonances.size() > 2)
363  {
364  for (unsigned int p = 3; p <= v_intermediateResonances.size(); ++p) dummyTrackList = findNProngs(v_intermediateResonances,
365  dummyTrackID, dummyTrackList,
366  (int) v_intermediateResonances.size(), (int) p);
367  }
368 
369  if (dummyTrackList.size() != 0)
370  {
371  std::vector<int> goodTrack{*i_it};
372  goodTracksThatMeetIntermediates.push_back(goodTrack);
373  }
374  }
375  }
376  else
377  {
378  goodTracksThatMeet = findTwoProngs(daughterParticles, goodTrackIndex, num_remaining_tracks);
379 
380  for (int p = 3; p <= num_remaining_tracks; ++p)
381  {
382  goodTracksThatMeet = findNProngs(daughterParticles, goodTrackIndex, goodTracksThatMeet, num_remaining_tracks, p);
383  }
384 
385  for (unsigned int i = 0; i < goodTracksThatMeet.size(); ++i)
386  {
387  std::vector<KFParticle> v_intermediateResonances(intermediateResonances, intermediateResonances + m_num_intermediate_states);
388  std::vector<std::vector<int>> dummyTrackList;
389  std::vector<int> dummyTrackID; //I already have the track ids stored in goodTracksThatMeet[i]
390  for (unsigned int j = 0; j < goodTracksThatMeet[i].size(); ++j)
391  {
392  v_intermediateResonances.push_back(daughterParticles[goodTracksThatMeet[i][j]]);
393  }
394  for (unsigned int k = 0; k < v_intermediateResonances.size(); ++k)
395  {
396  dummyTrackID.push_back(k);
397  }
398  dummyTrackList = findTwoProngs(v_intermediateResonances, dummyTrackID, (int) v_intermediateResonances.size());
399  for (unsigned int p = 3; p <= v_intermediateResonances.size(); ++p)
400  {
401  dummyTrackList = findNProngs(v_intermediateResonances, dummyTrackID, dummyTrackList, (int) v_intermediateResonances.size(), (int) p);
402  }
403 
404  if (dummyTrackList.size() != 0) goodTracksThatMeetIntermediates.push_back(goodTracksThatMeet[i]);
405  }
406  }
407 
408  return goodTracksThatMeetIntermediates;
409 }
410 
411 float KFParticle_Tools::eventDIRA(KFParticle particle, KFParticle vertex)
412 {
413  TMatrixD flightVector(3, 1);
414  TMatrixD momVector(3, 1);
415  flightVector(0, 0) = particle.GetX() - vertex.GetX();
416  flightVector(1, 0) = particle.GetY() - vertex.GetY();
417  flightVector(2, 0) = particle.GetZ() - vertex.GetZ();
418 
419  momVector(0, 0) = particle.GetPx();
420  momVector(1, 0) = particle.GetPy();
421  momVector(2, 0) = particle.GetPz();
422 
423  TMatrixD momDotFD(1, 1); //Calculate momentum dot flight distance
424  momDotFD = TMatrixD(momVector, TMatrixD::kTransposeMult, flightVector);
425  float f_momDotFD = momDotFD(0, 0);
426 
427  TMatrixD sizeOfMom(1, 1); //Calculates the size of the momentum vector
428  sizeOfMom = TMatrixD(momVector, TMatrixD::kTransposeMult, momVector);
429  float f_sizeOfMom = sqrt(sizeOfMom(0, 0));
430 
431  TMatrixD sizeOfFD(1, 1); //Calculates the size of the flight distance vector
432  sizeOfFD = TMatrixD(flightVector, TMatrixD::kTransposeMult, flightVector);
433  float f_sizeOfFD = sqrt(sizeOfFD(0, 0));
434 
435  return f_momDotFD / (f_sizeOfMom * f_sizeOfFD);
436 }
437 
438 float KFParticle_Tools::flightDistanceChi2(KFParticle particle, KFParticle vertex)
439 {
440  TMatrixD flightVector(3, 1);
441  TMatrixD flightDistanceCovariance(3, 3);
442 
443  KFParticle kfp_vertex(vertex);
444 
445  flightVector(0, 0) = particle.GetX() - kfp_vertex.GetX();
446  flightVector(1, 0) = particle.GetY() - kfp_vertex.GetY();
447  flightVector(2, 0) = particle.GetZ() - kfp_vertex.GetZ();
448 
449  for (int i = 0; i < 3; i++)
450  {
451  for (int j = 0; j < 3; j++)
452  {
453  flightDistanceCovariance(i, j) = particle.GetCovariance(i, j) + kfp_vertex.GetCovariance(i, j);
454  }
455  }
456 
457  TMatrixD anInverseMatrix(3, 3);
458  anInverseMatrix = flightDistanceCovariance.Invert();
459  TMatrixD m_chi2Value(1, 1);
460  m_chi2Value = TMatrixD(flightVector, TMatrixD::kTransposeMult, anInverseMatrix * flightVector);
461  return m_chi2Value(0, 0);
462 }
463 
464 std::tuple<KFParticle, bool> KFParticle_Tools::buildMother(KFParticle vDaughters[], std::string daughterOrder[],
465  bool isIntermediate, int intermediateNumber, int nTracks,
466  bool constrainMass, float required_vertexID)
467 {
468  KFParticle mother;
469  KFParticle inputTracks[nTracks];
470 
471  mother.SetConstructMethod(2);
472 
473  bool daughterMassCheck = true;
474  float unique_vertexID = 0;
475  float daughterMass = 0;
476 
477  //Figure out if the decay has reco. tracks mixed with resonances
478  int num_tracks_used_by_intermediates = 0;
479  for (int i = 0; i < m_num_intermediate_states; ++i) num_tracks_used_by_intermediates += m_num_tracks_from_intermediate[i];
480  int num_remaining_tracks = m_num_tracks - num_tracks_used_by_intermediates;
481 
482  for (int i = 0; i < nTracks; ++i)
483  {
484  daughterMass = constrainMass ? particleMasses.find(daughterOrder[i].c_str())->second.second : vDaughters[i].GetMass();
485  if ((num_remaining_tracks > 0 && i >= m_num_intermediate_states) || isIntermediate)
486  {
487  daughterMass = particleMasses.find(daughterOrder[i].c_str())->second.second;
488  }
489  inputTracks[i].Create(vDaughters[i].Parameters(),
490  vDaughters[i].CovarianceMatrix(),
491  (Int_t) vDaughters[i].GetQ(),
492  daughterMass);
493  mother.AddDaughter(inputTracks[i]);
494  unique_vertexID += vDaughters[i].GetQ() * particleMasses.find(daughterOrder[i].c_str())->second.second;
495  }
496 
497  if (isIntermediate) mother.SetPDG(particleMasses.find(m_intermediate_name[intermediateNumber].c_str())->second.first);
498  if (!isIntermediate && !m_mother_name_Tools.empty()) mother.SetPDG(particleMasses.find(m_mother_name_Tools.c_str())->second.first);
499 
500  bool chargeCheck;
502  chargeCheck = std::abs(unique_vertexID) == std::abs(required_vertexID) ? 1 : 0;
503  else
504  chargeCheck = unique_vertexID == required_vertexID ? 1 : 0;
505 
506  for (int j = 0; j < nTracks; ++j)
507  {
508  inputTracks[j].SetProductionVertex(mother);
510  {
511  if (inputTracks[j].GetMass() == 0) daughterMassCheck = false;
512  }
513  }
514 
515  float calculated_mass, calculated_mass_err;
516  mother.GetMass(calculated_mass, calculated_mass_err);
517  float calculated_pt = mother.GetPt();
518 
519  float min_mass = isIntermediate ? m_intermediate_mass_range[intermediateNumber].first : m_min_mass;
520  float max_mass = isIntermediate ? m_intermediate_mass_range[intermediateNumber].second : m_max_mass;
521  float min_pt = isIntermediate ? m_intermediate_min_pt[intermediateNumber] : m_mother_pt;
522 
523  bool goodCandidate = false;
524  if (calculated_mass >= min_mass && calculated_mass <= max_mass &&
525  calculated_pt >= min_pt && daughterMassCheck && chargeCheck)
526  goodCandidate = true;
527 
528  // Check the requirements of an intermediate states against this mother and re-do goodCandidate
529  if (goodCandidate && m_has_intermediates && !isIntermediate) //The decay has intermediate states and we are now looking at the mother
530  {
531  for (int k = 0; k < m_num_intermediate_states; ++k)
532  {
533  float intermediate_DIRA = eventDIRA(vDaughters[k], mother);
534  float intermediate_FDchi2 = flightDistanceChi2(vDaughters[k], mother);
535  if (intermediate_DIRA < m_intermediate_min_dira[k] ||
536  intermediate_FDchi2 < m_intermediate_min_fdchi2[k])
537  goodCandidate = false;
538  }
539  }
540 
541  return std::make_tuple(mother, goodCandidate);
542 }
543 
544 void KFParticle_Tools::constrainToVertex(KFParticle &particle, bool &goodCandidate, KFParticle &vertex)
545 {
546  KFParticle particleCopy = particle;
547  particleCopy.SetProductionVertex(vertex);
548  particleCopy.TransportToDecayVertex();
549 
550  float calculated_decayTime;
551  float calculated_decayTimeErr;
552  float calculated_decayLength;
553  float calculated_decayLengthErr;
554 
555  particleCopy.GetLifeTime(calculated_decayTime, calculated_decayTimeErr);
556  particleCopy.GetDecayLength(calculated_decayLength, calculated_decayLengthErr);
557 
558  float calculated_fdchi2 = flightDistanceChi2(particle, vertex);
559  float calculated_dira = eventDIRA(particle, vertex);
560  float calculated_ipchi2 = particle.GetDeviationFromVertex(vertex);
561 
562  goodCandidate = false;
563 
564  const float speed = 2.99792458e-1;
565  calculated_decayTime /= speed;
566 
567  if (calculated_fdchi2 >= m_fdchi2 && calculated_ipchi2 <= m_mother_ipchi2
568  && isInRange(m_dira_min, calculated_dira, m_dira_max)
569  && isInRange(m_min_decayTime, calculated_decayTime, m_max_decayTime)
570  && isInRange(m_min_decayLength, calculated_decayLength, m_max_decayLength))
571  goodCandidate = true;
572 }
573 
574 std::tuple<KFParticle, bool> KFParticle_Tools::getCombination(KFParticle vDaughters[], std::string daughterOrder[], KFParticle vertex, bool constrain_to_vertex, bool isIntermediate, int intermediateNumber, int nTracks, bool constrainMass, float required_vertexID)
575 {
576  KFParticle candidate;
577  bool isGoodCandidate;
578 
579  std::tie(candidate, isGoodCandidate) = buildMother(vDaughters, daughterOrder, isIntermediate, intermediateNumber, nTracks, constrainMass, required_vertexID);
580 
581  if (constrain_to_vertex && isGoodCandidate && !isIntermediate) constrainToVertex(candidate, isGoodCandidate, vertex);
582 
583  return std::make_tuple(candidate, isGoodCandidate);
584 }
585 
586 std::vector<std::vector<std::string>> KFParticle_Tools::findUniqueDaughterCombinations(int start, int end)
587 {
588  std::vector<int> vect_permutations;
589  std::vector<std::vector<std::string>> uniqueCombinations;
590  std::map<int, std::string> daughterMap;
591  for (int i = start; i < end; i++)
592  {
593  daughterMap.insert(std::pair<int, std::string>(i, m_daughter_name[i]));
594  vect_permutations.push_back(i);
595  }
596  int *permutations = &vect_permutations[0];
597 
598  do
599  {
600  std::vector<std::string> combination;
601  for (int i = 0; i < (end - start); i++) combination.push_back(daughterMap.find(permutations[i])->second);
602  uniqueCombinations.push_back(combination);
603  } while (std::next_permutation(permutations, permutations + vect_permutations.size()));
604 
605  removeDuplicates(uniqueCombinations);
606 
607  return uniqueCombinations;
608 }
609 
610 double KFParticle_Tools::calculateEllipsoidRadius(int posOrNeg, double sigma_ii, double sigma_jj, double sigma_ij)
611 { //Note - Only works for a 2D ellipsoid OR rotated nD ellipsoid to avoid projections
612  if (std::abs(posOrNeg) != 1)
613  {
614  std::cout << "You have set posOrNeg to " << posOrNeg << ". This value must be +/- 1! Skipping" << std::endl;
615  return 0;
616  }
617 
618  double r_ij = sqrt((sigma_ii + sigma_jj) / 2 + posOrNeg * (sqrt(pow((sigma_ii - sigma_jj) / 2, 2) + pow(sigma_ij, 2))));
619 
620  return r_ij;
621 }
622 
624 {
625  TMatrixD cov_matrix(3, 3);
626 
627  for (int i = 0; i < 3; ++i)
628  for (int j = 0; j < 3; ++j)
629  cov_matrix(i, j) = particle.GetCovariance(i, j);
630 
631  float volume;
632  if (cov_matrix(0, 0) * cov_matrix(1, 1) * cov_matrix(2, 2) == 0)
633  volume = 0;
634  else
635  volume = (4. / 3.) * M_PI * sqrt((std::abs(cov_matrix.Determinant()))); //The covariance matrix is error-squared
636 
637  return volume;
638 }
639 
640 float KFParticle_Tools::calculateJT(KFParticle mother, KFParticle daughter)
641 {
642  Eigen::Vector3f motherP = Eigen::Vector3f(mother.GetPx(), mother.GetPy(), mother.GetPz());
643  Eigen::Vector3f daughterP = Eigen::Vector3f(daughter.GetPx(), daughter.GetPy(), daughter.GetPz());
644 
645  Eigen::Vector3f motherP_X_daughterP = motherP.cross(daughterP);
646 
647  float jT = (motherP_X_daughterP.norm()) / motherP.norm();
648 
649  return jT;
650 }
651 
652 bool KFParticle_Tools::isInRange(float min, float value, float max)
653 {
654  return min <= value && value <= max;
655 }
656 
657 void KFParticle_Tools::removeDuplicates(std::vector<double> &v)
658 {
659  auto end = v.end();
660  for (auto it = v.begin(); it != end; ++it)
661  {
662  end = remove(it + 1, end, *it);
663  }
664  v.erase(end, v.end());
665 }
666 
667 void KFParticle_Tools::removeDuplicates(std::vector<int> &v)
668 {
669  auto end = v.end();
670  for (auto it = v.begin(); it != end; ++it)
671  {
672  end = remove(it + 1, end, *it);
673  }
674  v.erase(end, v.end());
675 }
676 
677 void KFParticle_Tools::removeDuplicates(std::vector<std::vector<int>> &v)
678 {
679  auto end = v.end();
680  for (auto it = v.begin(); it != end; ++it)
681  {
682  end = remove(it + 1, end, *it);
683  }
684  v.erase(end, v.end());
685 }
686 
687 void KFParticle_Tools::removeDuplicates(std::vector<std::vector<std::string>> &v)
688 {
689  auto end = v.end();
690  for (auto it = v.begin(); it != end; ++it)
691  {
692  end = remove(it + 1, end, *it);
693  }
694  v.erase(end, v.end());
695 }
696 
698 {
699  std::cout << "Track ID: " << particle.Id() << std::endl;
700  std::cout << "PDG ID: " << particle.GetPDG() << ", charge: " << (int) particle.GetQ() << ", mass: " << particle.GetMass() << " GeV" << std::endl;
701  std::cout << "(px,py,pz) = (" << particle.GetPx() << " +/- " << sqrt(particle.GetCovariance(3, 3)) << ", ";
702  std::cout << particle.GetPy() << " +/- " << sqrt(particle.GetCovariance(4, 4)) << ", ";
703  std::cout << particle.GetPz() << " +/- " << sqrt(particle.GetCovariance(5, 5)) << ") GeV" << std::endl;
704  std::cout << "(x,y,z) = (" << particle.GetX() << " +/- " << sqrt(particle.GetCovariance(0, 0)) << ", ";
705  std::cout << particle.GetY() << " +/- " << sqrt(particle.GetCovariance(1, 1)) << ", ";
706  std::cout << particle.GetZ() << " +/- " << sqrt(particle.GetCovariance(2, 2)) << ") cm\n"
707  << std::endl;
708 }