11 #include <tracking/trackFindingCDC/mclookup/CDCMCTrackStore.h>
13 #include <tracking/trackFindingCDC/mclookup/CDCSimHitLookUp.h>
14 #include <tracking/trackFindingCDC/mclookup/CDCMCMap.h>
15 #include <tracking/trackFindingCDC/mclookup/CDCMCManager.h>
17 #include <tracking/trackFindingCDC/topology/CDCWireTopology.h>
19 #include <tracking/trackFindingCDC/utilities/Functional.h>
20 #include <tracking/trackFindingCDC/utilities/Algorithms.h>
22 #include <cdc/dataobjects/CDCHit.h>
23 #include <cdc/dataobjects/CDCSimHit.h>
24 #include <mdst/dataobjects/MCParticle.h>
27 using namespace TrackFindingCDC;
38 B2DEBUG(200,
"In CDCMCTrackStore::clear()");
57 B2DEBUG(200,
"In CDCMCTrackStore::fill()");
81 B2DEBUG(100,
"m_inTrackIds.size(): " <<
m_inTrackIds.size());
84 B2DEBUG(100,
"m_nLoops.size(): " <<
m_nLoops.size());
92 B2WARNING(
"CDCMCMap not set. Cannot create tracks");
100 const MCParticle* ptrMCParticle = std::get<const MCParticle* const>(relation);
101 const CDCHit* ptrHit = std::get<const CDCHit*>(relation);
105 ITrackType mcParticleIdx = ptrMCParticle->getArrayIndex();
127 ITrackType mcParticleIdx = mcTrackAndMCParticleIdx.first;
130 if (mcTrack.empty())
continue;
135 auto superLayerRanges = adjacent_groupby(mcTrack.begin(), mcTrack.end(), [](
const CDCHit * hit) {
136 return hit->getISuperLayer();
139 std::vector<CDCHitVector> superLayerRuns;
140 for (
const auto& superLayerRange : superLayerRanges) {
141 superLayerRuns.push_back({superLayerRange.begin(), superLayerRange.end()});
144 std::vector<std::vector<CDCHitVector>::iterator> smallSuperLayerRuns;
145 for (
auto itSuperLayerRun = superLayerRuns.begin();
146 itSuperLayerRun != superLayerRuns.end();
148 if (itSuperLayerRun->size() < 3) smallSuperLayerRuns.push_back(itSuperLayerRun);
152 for (
auto itSuperLayerRun : smallSuperLayerRuns) {
153 ISuperLayer iSL = itSuperLayerRun->front()->getISuperLayer();
156 auto itSuperLayerRunBefore = superLayerRuns.end();
157 int hitDistanceBefore = INT_MAX;
158 if (std::distance(superLayerRuns.begin(), itSuperLayerRun) >= 2) {
159 itSuperLayerRunBefore = itSuperLayerRun - 2;
160 if (itSuperLayerRunBefore->front()->getISuperLayer() == iSL) {
161 hitDistanceBefore = (itSuperLayerRunBefore - 1)->size();
163 itSuperLayerRunBefore = superLayerRuns.end();
167 auto itSuperLayerRunAfter = superLayerRuns.end();
168 int hitDistanceAfter = INT_MAX;
169 if (std::distance(itSuperLayerRun, superLayerRuns.end()) > 2) {
170 itSuperLayerRunAfter = itSuperLayerRun + 2;
171 if (itSuperLayerRunAfter->front()->getISuperLayer() == iSL) {
172 hitDistanceAfter = (itSuperLayerRunAfter + 1)->size();
174 itSuperLayerRunAfter = superLayerRuns.end();
178 auto itMergeSuperLayerRun = superLayerRuns.end();
179 bool mergeBefore =
false;
180 if (hitDistanceBefore < hitDistanceAfter) {
181 itMergeSuperLayerRun = itSuperLayerRunBefore;
184 itMergeSuperLayerRun = itSuperLayerRunAfter;
188 if (itMergeSuperLayerRun == superLayerRuns.end())
continue;
189 else if (mergeBefore) {
190 itMergeSuperLayerRun->insert(itMergeSuperLayerRun->end(), itSuperLayerRun->begin(), itSuperLayerRun->end());
191 itSuperLayerRun->clear();
193 itMergeSuperLayerRun->insert(itMergeSuperLayerRun->begin(), itSuperLayerRun->begin(), itSuperLayerRun->end());
194 itSuperLayerRun->clear();
199 erase_remove_if(superLayerRuns,
Size() == 0u);
203 if (lhs.empty() or rhs.empty())
return true;
204 if (lhs.front()->getISuperLayer() != rhs.front()->getISuperLayer())
return false;
205 lhs.insert(lhs.end(), rhs.begin(), rhs.end());
209 erase_unique(superLayerRuns, mergeSameSuperLayer);
212 for (
const CDCHitVector& superLayerRun : superLayerRuns) {
214 auto areNeighbors = [&wireTopology](
const CDCHit * lhs,
const CDCHit * rhs) {
216 WireID rhsWireID(rhs->getISuperLayer(), rhs->getILayer(), rhs->getIWire());
220 lhsWireID == rhsWireID);
223 auto segmentRanges = unique_ranges(superLayerRun.begin(), superLayerRun.end(), areNeighbors);
226 mcSegments.emplace_back(segmentRange.begin(), segmentRange.end());
240 B2WARNING(
"CDCMCMap not set. Cannot sort track");
246 std::stable_sort(mcTrack.begin(), mcTrack.end(),
247 [&simHitLookUp](
const CDCHit * ptrHit,
const CDCHit * ptrOtherHit) ->
bool {
249 const CDCSimHit* ptrSimHit = simHitLookUp.getClosestPrimarySimHit(ptrHit);
250 const CDCSimHit* ptrOtherSimHit = simHitLookUp.getClosestPrimarySimHit(ptrOtherHit);
254 B2FATAL(
"No CDCSimHit for CDCHit");
257 if (not ptrOtherSimHit)
259 B2FATAL(
"No CDCSimHit for CDCHit");
262 double secondaryFlightTime = ptrSimHit->getFlightTime();
263 double otherSecondaryFlightTime = ptrOtherSimHit->getFlightTime();
266 return (secondaryFlightTime < std::fmin(INFINITY, otherSecondaryFlightTime));
276 void CDCMCTrackStore::fillInTrackId()
279 for (
const std::pair<ITrackType, CDCHitVector>& mcTrackAndMCParticleIdx : getMCTracksByMCParticleIdx()) {
281 const CDCHitVector& mcTrack = mcTrackAndMCParticleIdx.second;
285 for (
const CDCHit* ptrHit : mcTrack) {
287 m_inTrackIds[ptrHit] = iHit;
293 void CDCMCTrackStore::fillInTrackSegmentId()
295 for (
const std::pair<ITrackType, std::vector<CDCHitVector> >& mcSegmentsAndMCParticleIdx : getMCSegmentsByMCParticleIdx()) {
296 const std::vector<CDCHitVector>& mcSegments = mcSegmentsAndMCParticleIdx.second;
301 for (
const CDCHit* ptrHit : mcSegment) {
303 m_inTrackSegmentIds[ptrHit] = iSegment;
310 void CDCMCTrackStore::fillNLoopsAndNPassedSuperLayers()
313 for (
const std::pair<ITrackType, std::vector<CDCHitVector> >& mcSegmentsAndMCParticleIdx : getMCSegmentsByMCParticleIdx()) {
314 const std::vector<CDCHitVector>& mcSegments = mcSegmentsAndMCParticleIdx.second;
317 int nPassedSuperLayers = 0;
321 if (ptrLastMCSegment and changedSuperLayer(*ptrLastMCSegment, mcSegment)) {
322 ++nPassedSuperLayers;
326 if (ptrLastMCSegment->front()->getISuperLayer() == 0 and
327 mcSegment.front()->getISuperLayer() == 0) {
332 for (
const CDCHit* ptrHit : mcSegment) {
333 m_nPassedSuperLayers[ptrHit] = nPassedSuperLayers;
334 m_nLoops[ptrHit] = nLoops;
337 ptrLastMCSegment = &mcSegment;
346 const CDCHit* ptrHit = mcSegment.front();
347 const CDCHit* ptrNextHit = nextMCSegment.front();
352 const CDCHit& hit = *ptrHit;
353 const CDCHit& nextHit = *ptrNextHit;
357 }
else if (hit.getISuperLayer() == 0) {
366 if (pos.dotXY(nextPos) < 0)
return true;
367 if ((nextPos - pos).dotXY(nextMom) < 0)
return true;
368 if ((nextPos - pos).dotXY(mom) < 0)
return true;
378 Index CDCMCTrackStore::getInTrackId(
const CDCHit* ptrHit)
const
381 auto itFoundHit = m_inTrackIds.find(ptrHit);
382 return itFoundHit == m_inTrackIds.end() ? c_InvalidIndex : itFoundHit->second;
388 Index CDCMCTrackStore::getInTrackSegmentId(
const CDCHit* ptrHit)
const
391 auto itFoundHit = m_inTrackSegmentIds.find(ptrHit);
392 return itFoundHit == m_inTrackSegmentIds.end() ? c_InvalidIndex : itFoundHit->second;
397 Index CDCMCTrackStore::getNPassedSuperLayers(
const CDCHit* ptrHit)
const
400 auto itFoundHit = m_nPassedSuperLayers.find(ptrHit);
401 return itFoundHit == m_nPassedSuperLayers.end() ? c_InvalidIndex : itFoundHit->second;
405 Index CDCMCTrackStore::getNLoops(
const CDCHit* ptrHit)
const
408 auto itFoundHit = m_nLoops.find(ptrHit);
409 return itFoundHit == m_nLoops.end() ? c_InvalidIndex : itFoundHit->second;