38 #include <KFParticle.h>
39 #include <KFParticleDatabase.h>
45 #include <Eigen/Dense>
56 : m_has_intermediates(
false)
67 , m_track_chi2ndof(4.)
69 , m_vertex_chi2ndof(15.)
75 , m_get_charge_conjugate(
true)
76 , m_vtx_map_node_name(
"SvtxVertexMap")
77 , m_trk_map_node_name(
"SvtxTrackMap")
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)
100 KFParticle kfp_vertex;
101 kfp_vertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
111 if (vertexMapName.empty())
114 vtxMN = vertexMapName;
116 std::vector<KFParticle> primaryVertices;
118 unsigned int vertexID = 0;
124 primaryVertices.push_back(
makeVertex(topNode));
125 primaryVertices[vertexID].SetId(iter->first);
129 return primaryVertices;
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)
150 KFParticle kfp_particle;
161 std::vector<KFParticle> daughterParticles;
163 unsigned int trackID = 0;
169 daughterParticles[trackID].SetId(iter->first);
173 return daughterParticles;
179 if (vertexMapName.empty())
182 vtxMN = vertexMapName;
189 return associatedVertex->size_tracks();
194 bool goodTrack =
false;
197 float min_ipchi2 = 0;
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);
214 float& minimumIP,
float& minimumIPchi2)
216 std::vector<float> ip, ipchi2;
218 for (
unsigned int i_verts = 0; i_verts < PVs.size(); ++i_verts)
220 ip.push_back(track.GetDistanceFromVertex(PVs[i_verts]));
221 ipchi2.push_back(track.GetDeviationFromVertex(PVs[i_verts]));
224 auto minmax_ip = minmax_element(ip.begin(), ip.end());
225 minimumIP = *minmax_ip.first;
226 auto minmax_ipchi2 = minmax_element(ipchi2.begin(), ipchi2.end());
227 minimumIPchi2 = *minmax_ipchi2.first;
234 std::vector<int> goodTrackIndex;
236 for (
unsigned int i_parts = 0; i_parts < daughterParticles.size(); ++i_parts)
238 if (
isGoodTrack(daughterParticles[i_parts], primaryVertices)) goodTrackIndex.push_back(i_parts);
243 return goodTrackIndex;
248 std::vector<std::vector<int>> goodTracksThatMeet;
250 for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
252 for (std::vector<int>::iterator j_it = goodTrackIndex.begin(); j_it != goodTrackIndex.end(); ++j_it)
256 if (daughterParticles[*i_it].GetDistanceFromParticle(daughterParticles[*j_it]) <=
m_comb_DCA)
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};
266 goodTracksThatMeet.push_back(combination);
274 goodTracksThatMeet.push_back(combination);
281 return goodTracksThatMeet;
285 std::vector<int> goodTrackIndex,
286 std::vector<std::vector<int>> goodTracksThatMeet,
287 int nRequiredTracks,
unsigned int nProngs)
289 unsigned int nGoodProngs = goodTracksThatMeet.size();
291 for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
293 for (
unsigned int i_prongs = 0; i_prongs < nGoodProngs; ++i_prongs)
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)
301 for (
unsigned int i = 0; i < nProngs - 1; ++i)
303 if (daughterParticles[*i_it].GetDistanceFromParticle(daughterParticles[goodTracksThatMeet[i_prongs][i]]) >
m_comb_DCA)
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)
317 particleVertex += daughterParticles[goodTracksThatMeet[i_prongs][i]];
318 combination.push_back(goodTracksThatMeet[i_prongs][i]);
320 float vertexchi2ndof = particleVertex.GetChi2() / particleVertex.GetNDF();
322 if ((
unsigned int) nRequiredTracks == nProngs && vertexchi2ndof <=
m_vertex_chi2ndof)
324 goodTracksThatMeet.push_back(combination);
326 else if ((
unsigned int) nRequiredTracks == nProngs && vertexchi2ndof >
m_vertex_chi2ndof)
332 goodTracksThatMeet.push_back(combination);
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());
343 return goodTracksThatMeet;
348 std::vector<std::vector<int>> goodTracksThatMeet, goodTracksThatMeetIntermediates;
349 if (num_remaining_tracks == 1)
351 for (std::vector<int>::iterator i_it = goodTrackIndex.begin(); i_it != goodTrackIndex.end(); ++i_it)
353 std::vector<KFParticle> v_intermediateResonances(intermediateResonances, intermediateResonances +
m_num_intermediate_states);
354 std::vector<std::vector<int>> dummyTrackList;
355 std::vector<int> dummyTrackID;
356 v_intermediateResonances.insert(end(v_intermediateResonances), daughterParticles[*i_it]);
357 for (
unsigned int k = 0;
k < v_intermediateResonances.size(); ++
k)
359 dummyTrackID.push_back(
k);
361 dummyTrackList =
findTwoProngs(v_intermediateResonances, dummyTrackID, (
int) v_intermediateResonances.size());
362 if (v_intermediateResonances.size() > 2)
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);
369 if (dummyTrackList.size() != 0)
371 std::vector<int> goodTrack{*i_it};
372 goodTracksThatMeetIntermediates.push_back(goodTrack);
378 goodTracksThatMeet =
findTwoProngs(daughterParticles, goodTrackIndex, num_remaining_tracks);
380 for (
int p = 3;
p <= num_remaining_tracks; ++
p)
382 goodTracksThatMeet =
findNProngs(daughterParticles, goodTrackIndex, goodTracksThatMeet, num_remaining_tracks,
p);
385 for (
unsigned int i = 0; i < goodTracksThatMeet.size(); ++i)
387 std::vector<KFParticle> v_intermediateResonances(intermediateResonances, intermediateResonances +
m_num_intermediate_states);
388 std::vector<std::vector<int>> dummyTrackList;
389 std::vector<int> dummyTrackID;
390 for (
unsigned int j = 0; j < goodTracksThatMeet[i].size(); ++j)
392 v_intermediateResonances.push_back(daughterParticles[goodTracksThatMeet[i][j]]);
394 for (
unsigned int k = 0;
k < v_intermediateResonances.size(); ++
k)
396 dummyTrackID.push_back(
k);
398 dummyTrackList =
findTwoProngs(v_intermediateResonances, dummyTrackID, (
int) v_intermediateResonances.size());
399 for (
unsigned int p = 3;
p <= v_intermediateResonances.size(); ++
p)
401 dummyTrackList =
findNProngs(v_intermediateResonances, dummyTrackID, dummyTrackList, (
int) v_intermediateResonances.size(), (
int)
p);
404 if (dummyTrackList.size() != 0) goodTracksThatMeetIntermediates.push_back(goodTracksThatMeet[i]);
408 return goodTracksThatMeetIntermediates;
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();
419 momVector(0, 0) = particle.GetPx();
420 momVector(1, 0) = particle.GetPy();
421 momVector(2, 0) = particle.GetPz();
423 TMatrixD momDotFD(1, 1);
424 momDotFD = TMatrixD(momVector, TMatrixD::kTransposeMult, flightVector);
425 float f_momDotFD = momDotFD(0, 0);
427 TMatrixD sizeOfMom(1, 1);
428 sizeOfMom = TMatrixD(momVector, TMatrixD::kTransposeMult, momVector);
429 float f_sizeOfMom = sqrt(sizeOfMom(0, 0));
431 TMatrixD sizeOfFD(1, 1);
432 sizeOfFD = TMatrixD(flightVector, TMatrixD::kTransposeMult, flightVector);
433 float f_sizeOfFD = sqrt(sizeOfFD(0, 0));
435 return f_momDotFD / (f_sizeOfMom * f_sizeOfFD);
440 TMatrixD flightVector(3, 1);
441 TMatrixD flightDistanceCovariance(3, 3);
443 KFParticle kfp_vertex(vertex);
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();
449 for (
int i = 0; i < 3; i++)
451 for (
int j = 0; j < 3; j++)
453 flightDistanceCovariance(i, j) = particle.GetCovariance(i, j) + kfp_vertex.GetCovariance(i, j);
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);
465 bool isIntermediate,
int intermediateNumber,
int nTracks,
466 bool constrainMass,
float required_vertexID)
469 KFParticle inputTracks[nTracks];
471 mother.SetConstructMethod(2);
473 bool daughterMassCheck =
true;
474 float unique_vertexID = 0;
475 float daughterMass = 0;
478 int num_tracks_used_by_intermediates = 0;
480 int num_remaining_tracks =
m_num_tracks - num_tracks_used_by_intermediates;
482 for (
int i = 0; i < nTracks; ++i)
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)
487 daughterMass = particleMasses.find(daughterOrder[i].c_str())->second.second;
489 inputTracks[i].Create(vDaughters[i].
Parameters(),
490 vDaughters[i].CovarianceMatrix(),
491 (
Int_t) vDaughters[i].GetQ(),
493 mother.AddDaughter(inputTracks[i]);
494 unique_vertexID += vDaughters[i].GetQ() * particleMasses.find(daughterOrder[i].c_str())->second.second;
502 chargeCheck =
std::abs(unique_vertexID) ==
std::abs(required_vertexID) ? 1 : 0;
504 chargeCheck = unique_vertexID == required_vertexID ? 1 : 0;
506 for (
int j = 0; j < nTracks; ++j)
508 inputTracks[j].SetProductionVertex(mother);
511 if (inputTracks[j].GetMass() == 0) daughterMassCheck =
false;
515 float calculated_mass, calculated_mass_err;
516 mother.GetMass(calculated_mass, calculated_mass_err);
517 float calculated_pt = mother.GetPt();
523 bool goodCandidate =
false;
524 if (calculated_mass >= min_mass && calculated_mass <= max_mass &&
525 calculated_pt >= min_pt && daughterMassCheck && chargeCheck)
526 goodCandidate =
true;
533 float intermediate_DIRA =
eventDIRA(vDaughters[
k], mother);
537 goodCandidate =
false;
541 return std::make_tuple(mother, goodCandidate);
547 particleCopy.SetProductionVertex(vertex);
548 particleCopy.TransportToDecayVertex();
550 float calculated_decayTime;
551 float calculated_decayTimeErr;
552 float calculated_decayLength;
553 float calculated_decayLengthErr;
555 particleCopy.GetLifeTime(calculated_decayTime, calculated_decayTimeErr);
556 particleCopy.GetDecayLength(calculated_decayLength, calculated_decayLengthErr);
559 float calculated_dira =
eventDIRA(particle, vertex);
560 float calculated_ipchi2 = particle.GetDeviationFromVertex(vertex);
562 goodCandidate =
false;
564 const float speed = 2.99792458e-1;
565 calculated_decayTime /= speed;
571 goodCandidate =
true;
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)
576 KFParticle candidate;
577 bool isGoodCandidate;
579 std::tie(candidate, isGoodCandidate) =
buildMother(vDaughters, daughterOrder, isIntermediate, intermediateNumber, nTracks, constrainMass, required_vertexID);
581 if (constrain_to_vertex && isGoodCandidate && !isIntermediate)
constrainToVertex(candidate, isGoodCandidate, vertex);
583 return std::make_tuple(candidate, isGoodCandidate);
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++)
593 daughterMap.insert(std::pair<int, std::string>(i,
m_daughter_name[i]));
594 vect_permutations.push_back(i);
596 int *permutations = &vect_permutations[0];
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()));
607 return uniqueCombinations;
614 std::cout <<
"You have set posOrNeg to " << posOrNeg <<
". This value must be +/- 1! Skipping" << std::endl;
618 double r_ij = sqrt((sigma_ii + sigma_jj) / 2 + posOrNeg * (sqrt(pow((sigma_ii - sigma_jj) / 2, 2) + pow(sigma_ij, 2))));
625 TMatrixD cov_matrix(3, 3);
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);
632 if (cov_matrix(0, 0) * cov_matrix(1, 1) * cov_matrix(2, 2) == 0)
635 volume = (4. / 3.) *
M_PI * sqrt((
std::abs(cov_matrix.Determinant())));
642 Eigen::Vector3f motherP = Eigen::Vector3f(mother.GetPx(), mother.GetPy(), mother.GetPz());
643 Eigen::Vector3f daughterP = Eigen::Vector3f(daughter.GetPx(), daughter.GetPy(), daughter.GetPz());
645 Eigen::Vector3f motherP_X_daughterP = motherP.cross(daughterP);
647 float jT = (motherP_X_daughterP.norm()) / motherP.norm();
654 return min <= value && value <=
max;
660 for (
auto it = v.begin();
it != end; ++
it)
662 end =
remove(
it + 1, end, *
it);
664 v.erase(end, v.end());
670 for (
auto it = v.begin();
it != end; ++
it)
672 end =
remove(
it + 1, end, *
it);
674 v.erase(end, v.end());
680 for (
auto it =
v.begin();
it != end; ++
it)
682 end =
remove(
it + 1, end, *
it);
684 v.erase(end,
v.end());
690 for (
auto it =
v.begin();
it != end; ++
it)
692 end =
remove(
it + 1, end, *
it);
694 v.erase(end,
v.end());
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"