332 B2DEBUG(20,
"Skipping MC Track Finder as there are no MC Particles registered in the DataStore.");
337 const int eventCounter = eventMetaDataPtr->getEvent();
338 B2DEBUG(20,
"******* MCTrackFinderModule processing event number: " << eventCounter <<
" *******");
342 B2DEBUG(20,
"MCTrackFinder: total Number of MCParticles: " << nMcParticles);
346 B2DEBUG(20,
"MCTrackFinder: Number of PXDTrueHits: " << nPXDHits);
349 const int nMcPartToPXDHits = mcPartToPXDTrueHits.
getEntries();
350 B2DEBUG(20,
"MCTrackFinder: Number of relations between MCParticles and PXDHits: " << nMcPartToPXDHits);
354 B2DEBUG(20,
"MCTrackFinder: Number of PXDClusters: " << nPXDClusters);
357 const int nPxdClusterToMCPart = pxdClusterToMCParticle.
getEntries();
358 B2DEBUG(20,
"MCTrackFinder: Number of relations between PXDCluster and MCParticles: " << nPxdClusterToMCPart);
362 B2DEBUG(20,
"MCTrackFinder: Number of SVDDHits: " << nSVDHits);
365 const int nMcPartToSVDHits = mcPartToSVDTrueHits.
getEntries();
366 B2DEBUG(20,
"MCTrackFinder: Number of relations between MCParticles and SVDHits: " << nMcPartToSVDHits);
370 B2DEBUG(20,
"MCTrackFinder: Number of SVDClusters: " << nSVDClusters);
373 const int nSvdClusterToMCPart = svdClusterToMCParticle.
getEntries();
374 B2DEBUG(20,
"MCTrackFinder: Number of relations between SVDCluster and MCParticles: " << nSvdClusterToMCPart);
377 const int nCDCHits =
m_CDCHits.getEntries();
378 B2DEBUG(20,
"MCTrackFinder: Number of CDCHits: " << nCDCHits);
381 const int nMcPartToCDCHits = mcPartToCDCHits.
getEntries();
382 B2DEBUG(20,
"MCTrackFinder: Number of relations between MCParticles and CDCHits: " << nMcPartToCDCHits);
385 B2DEBUG(20,
"MCTrackFinder: Number of CDCSimHits: " << nCDCSimHits);
388 const int nCdcSimHitToHitRel = cdcSimHitToHitRel.
getEntries();
389 B2DEBUG(20,
"MCTrackFinder: Number of relations between CDCSimHit and CDCHits: " << nCdcSimHitToHitRel);
395 std::set<int> alreadyConsumedMCParticles;
396 for (
int iPart = 0; iPart < nMcParticles; ++iPart) {
397 if (alreadyConsumedMCParticles.count(iPart))
continue;
398 alreadyConsumedMCParticles.insert(iPart);
403 B2DEBUG(20,
"Particle that did not propagate significantly cannot make track.");
408 int mcParticleProperties = 0;
410 mcParticleProperties += 1;
413 mcParticleProperties += 2;
416 mcParticleProperties += 4;
419 mcParticleProperties += 8;
422 mcParticleProperties += 16;
425 mcParticleProperties += 32;
428 mcParticleProperties += 64;
431 mcParticleProperties += 128;
435 B2DEBUG(20,
"PDG: " << aMcParticlePtr->
getPDG() <<
" | property mask of particle " << mcParticleProperties <<
441 B2DEBUG(20,
"particle energy too low. MC particle will be skipped");
448 if (nPdgCodes not_eq 0) {
449 const int currentPdgCode = aMcParticlePtr->
getPDG();
450 int nFalsePdgCodes = 0;
451 for (
int i = 0; i not_eq nPdgCodes; ++i) {
456 if (nFalsePdgCodes == nPdgCodes) {
457 B2DEBUG(20,
"particle does not have one of the user provided pdg codes and will therefore be skipped");
465 if (nFromPdgCodes not_eq 0) {
467 int nFalsePdgCodes = 0;
469 bool foundParent =
false;
470 while (currentParent not_eq
nullptr) {
471 int currentParentPdgCode = currentParent->
getPDG();
472 for (
int i = 0; i not_eq nFromPdgCodes; ++i) {
484 currentParent =
nullptr;
486 currentParent = currentParent->
getMother();
491 B2DEBUG(20,
"particle does not have and ancestor with one of the user provided pdg codes and will therefore be skipped");
498 if (abs(aMcParticlePtr->
getPDG()) > 1000000000
500 B2DEBUG(20,
"Skipped Baryon.");
506 if (!
m_neutrals && (aMcParticlePtr->
getCharge() == 0 && abs(aMcParticlePtr->
getPDG()) != Monopoles::c_monopolePDGCode)) {
507 B2DEBUG(20,
"particle does not have the right charge. MC particle will be skipped");
512 B2DEBUG(20,
"Build a track for the MCParticle with index: " << iPart <<
" (PDG: " << aMcParticlePtr->
getPDG() <<
")");
514 std::vector<const MCParticle*> trackMCParticles;
515 trackMCParticles.push_back(aMcParticlePtr);
518 std::vector<MCParticle*> daughters = aMcParticlePtr->
getDaughters();
519 std::vector<MCParticle*> decayInFlightParticles;
521 if (daughter->getSecondaryPhysicsProcess() > 200 and daughter->getCharge() != 0) {
522 decayInFlightParticles.push_back(daughter);
525 if (decayInFlightParticles.size() == 1) {
526 trackMCParticles.push_back(decayInFlightParticles.front());
527 alreadyConsumedMCParticles.insert(trackMCParticles.back()->getArrayIndex());
536 typedef std::tuple<double, int, RecoHitInformation::OriginTrackFinder, Const::EDetector> TimeHitIDDetector;
537 std::vector<TimeHitIDDetector> hitsWithTimeAndDetectorInformation;
544 for (
const MCParticle* trackMCParticle : trackMCParticles) {
547 for (
size_t i = 0; i < relatedClusters.
size(); ++i) {
548 bool isReassigned = relatedClusters.
weight(i) < 0;
552 auto mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderPriorityHit;
556 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
561 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
566 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
570 if (
m_discardAuxiliaryHits and mcFinder == RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit)
continue;
582 for (
const PXDTrueHit& pxdTrueHit : relatedTrueHits) {
586 if (std::find_if(relatedMCParticles.
begin(),
587 relatedMCParticles.
end(),
588 [trackMCParticle](
const MCParticle & mcParticle) {
589 return &mcParticle == trackMCParticle;
590 }) != relatedMCParticles.
end()) {
591 time = pxdTrueHit.getGlobalTime();
595 if (not std::isnan(time)) {
596 hitsWithTimeAndDetectorInformation.emplace_back(time, pxdCluster->
getArrayIndex(), mcFinder, Const::PXD);
602 B2DEBUG(20,
" add " << hitCounter <<
" PXDClusters. " << relatedClusters.
size() - hitCounter <<
603 " PXDClusters were not added because they do not have a corresponding PXDTrueHit");
615 for (
const MCParticle* trackMCParticle : trackMCParticles) {
618 for (
size_t i = 0; i < relatedClusters.
size(); ++i) {
619 bool isReassigned = relatedClusters.
weight(i) < 0;
622 auto mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderPriorityHit;
626 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
631 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
636 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
640 if (
m_discardAuxiliaryHits and mcFinder == RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit)
continue;
652 for (
const SVDTrueHit& svdTrueHit : relatedTrueHits) {
656 if (std::find_if(relatedMCParticles.
begin(),
657 relatedMCParticles.
end(),
658 [trackMCParticle](
const MCParticle & mcParticle) {
659 return &mcParticle == trackMCParticle;
660 }) != relatedMCParticles.
end()) {
661 time = svdTrueHit.getGlobalTime();
665 if (not std::isnan(time)) {
666 hitsWithTimeAndDetectorInformation.emplace_back(time, svdCluster->
getArrayIndex(), mcFinder, Const::SVD);
672 B2DEBUG(20,
" add " << hitCounter <<
" SVDClusters. " << relatedClusters.
size() - hitCounter <<
673 " SVDClusters were not added because they do not have a corresponding SVDTrueHit");
686 std::array<int, 9> nHitsBySuperLayerId{};
688 for (
const MCParticle* trackMCParticle : trackMCParticles) {
691 for (
size_t i = 0; i < relatedHits.
size(); ++i) {
692 bool isReassigned = relatedHits.
weight(i) < 0;
699 unsigned short layerID = cdcHit->
getICLayer();
707 auto mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderPriorityHit;
711 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
716 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
721 mcFinder = RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit;
725 if (
m_discardAuxiliaryHits and mcFinder == RecoHitInformation::OriginTrackFinder::c_MCTrackFinderAuxiliaryHit)
continue;
729 if (not aCDCSimHitPtr) {
730 B2DEBUG(20,
" Skipping CDCHit without related CDCSimHit.");
737 if (superLayerId % 2 == 0) {
738 hitsWithTimeAndDetectorInformation.emplace_back(time, cdcHit->
getArrayIndex(), mcFinder, Const::CDC);
741 ++nHitsBySuperLayerId[superLayerId];
744 hitsWithTimeAndDetectorInformation.emplace_back(time, cdcHit->
getArrayIndex(), mcFinder, Const::CDC);
747 ++nHitsBySuperLayerId[superLayerId];
752 B2DEBUG(20,
" added " << nAxialHits <<
" axial and " << nStereoHits <<
" stereo CDCHits");
761 nHitsBySuperLayerId[0] == nAxialHits and
781 std::sort(hitsWithTimeAndDetectorInformation.begin(), hitsWithTimeAndDetectorInformation.end(),
782 [](
const TimeHitIDDetector & rhs,
const TimeHitIDDetector & lhs) {
783 return std::get<0>(rhs) < std::get<0>(lhs);
789 std::vector< std::vector<TimeHitIDDetector> > hitsWithTimeAndDetectorInformationVectors;
792 hitsWithTimeAndDetectorInformationVectors.push_back(hitsWithTimeAndDetectorInformation);
795 std::vector<TimeHitIDDetector>::size_type splitFromIdx = 0;
796 for (std::vector<TimeHitIDDetector>::size_type i = 1; i != hitsWithTimeAndDetectorInformation.size(); i++) {
798 double delta_t = (std::get<0>(hitsWithTimeAndDetectorInformation[i])
799 - std::get<0>(hitsWithTimeAndDetectorInformation[i - 1]));
803 hitsWithTimeAndDetectorInformationVectors
804 .emplace_back(hitsWithTimeAndDetectorInformation.begin() + splitFromIdx,
805 hitsWithTimeAndDetectorInformation.begin() + i);
810 hitsWithTimeAndDetectorInformationVectors
811 .emplace_back(hitsWithTimeAndDetectorInformation.begin() + splitFromIdx,
812 hitsWithTimeAndDetectorInformation.end());
817 B2DEBUG(20,
"We came pass all filter of the MCPartile and hit properties. TrackCandidate " << counter <<
818 " will be created from the MCParticle with index: " << iPart <<
" (PDG: " << aMcParticlePtr->
getPDG() <<
")");
824 ROOT::Math::XYZVector momentumTrue = aMcParticlePtr->
getMomentum();
828 ROOT::Math::XYZVector momentum = momentumTrue;
829 ROOT::Math::XYZVector position = positionTrue;
830 double time = timeTrue;
831 TVectorD stateSeed(6);
832 TMatrixDSym covSeed(6);
834 covSeed(0, 0) = 1; covSeed(1, 1) = 1; covSeed(2, 2) = 2 * 2;
835 covSeed(3, 3) = 0.1 * 0.1; covSeed(4, 4) = 0.1 * 0.1; covSeed(5, 5) = 0.2 * 0.2;
843 double smearedX = gRandom->Gaus(positionTrue.x(), smearing * positionTrue.x());
844 double smearedY = gRandom->Gaus(positionTrue.y(), smearing * positionTrue.y());
845 double smearedZ = gRandom->Gaus(positionTrue.z(), smearing * positionTrue.z());
846 position.SetXYZ(smearedX, smearedY, smearedZ);
847 double smearedPX = gRandom->Gaus(momentumTrue.x(), smearing * momentumTrue.x());
848 double smearedPY = gRandom->Gaus(momentumTrue.y(), smearing * momentumTrue.y());
849 double smearedPZ = gRandom->Gaus(momentumTrue.z(), smearing * momentumTrue.z());
850 momentum.SetXYZ(smearedPX, smearedPY, smearedPZ);
857 double smearedX = gRandom->Gaus(positionTrue.x(),
sqrt(
m_initialCov(0, 0)));
858 double smearedY = gRandom->Gaus(positionTrue.y(),
sqrt(
m_initialCov(1, 1)));
859 double smearedZ = gRandom->Gaus(positionTrue.z(),
sqrt(
m_initialCov(2, 2)));
860 position.SetXYZ(smearedX, smearedY, smearedZ);
861 double smearedPX = gRandom->Gaus(momentumTrue.x(),
sqrt(
m_initialCov(3, 3)));
862 double smearedPY = gRandom->Gaus(momentumTrue.y(),
sqrt(
m_initialCov(4, 4)));
863 double smearedPZ = gRandom->Gaus(momentumTrue.z(),
sqrt(
m_initialCov(5, 5)));
864 momentum.SetXYZ(smearedPX, smearedPY, smearedPZ);
870 for (
const auto& hitInformationVector : hitsWithTimeAndDetectorInformationVectors) {
873 short int charge =
static_cast<short int>(aMcParticlePtr->
getCharge());
876 if (hitInformationVector.size() != 0) {
878 time = std::get<0>(hitInformationVector.at(0));
881 const double beta_xy = momentum.Rho() / energy;
887 momentum = helix.getMomentumAtArcLength2D(arclength2D, Bz);
888 position = helix.getPositionAtArcLength2D(arclength2D);
900 B2DEBUG(20,
" --- Create relation between genfit::TrackCand " << counter <<
" and MCParticle " << iPart);
903 for (
const TimeHitIDDetector& hitInformation : hitInformationVector) {
906 const int hitID = std::get<1>(hitInformation);
907 const auto hitOriginMCFinderType = std::get<2>(hitInformation);
910 if (detectorInformation == Const::CDC) {
915 ROOT::Math::XYZVector simHitPos = aCDCSimHitPtr->
getPosTrack();
916 ROOT::Math::XYZVector simMom = aCDCSimHitPtr->
getMomentum();
917 ROOT::Math::XYZVector simHitPosOnWire = aCDCSimHitPtr->
getPosWire();
920 const unsigned short isRightHit = cdcGeometry.
getNewLeftRightRaw(simHitPosOnWire, simHitPos, simMom);
923 newRecoTrack->
addCDCHit(cdcHit, hitCounter, RecoHitInformation::RightLeftInformation::c_right, hitOriginMCFinderType);
925 newRecoTrack->
addCDCHit(cdcHit, hitCounter, RecoHitInformation::RightLeftInformation::c_left, hitOriginMCFinderType);
927 B2DEBUG(20,
"CDC hit " << hitID <<
" has reft/right sign " << isRightHit);
928 }
else if (detectorInformation == Const::PXD) {
930 newRecoTrack->
addPXDHit(pxdCluster, hitCounter, hitOriginMCFinderType);
931 }
else if (detectorInformation == Const::SVD) {
933 newRecoTrack->
addSVDHit(svdCluster, hitCounter, hitOriginMCFinderType);
938 B2DEBUG(20,
"New RecoTrack: #PXDHits: " << newRecoTrack->
getPXDHitList().size() <<