Belle II Software  release-05-01-25
CDCMCMap.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2014 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Oliver Frost *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/trackFindingCDC/mclookup/CDCMCMap.h>
11 
12 #include <framework/datastore/StoreArray.h>
13 
14 #include <cdc/dataobjects/CDCHit.h>
15 #include <cdc/dataobjects/CDCSimHit.h>
16 #include <mdst/dataobjects/MCParticle.h>
17 
18 #include <iostream>
19 
20 using namespace Belle2;
21 using namespace TrackFindingCDC;
22 
23 namespace {
24  template <class MapType>
25  void print_map(const MapType& m)
26  {
27  using const_iterator = typename MapType::const_iterator;
28  for (const_iterator iter = m.begin(), iend = m.end(); iter != iend; ++iter) {
29  std::cout << iter->first << "-->" << iter->second << std::endl;
30  }
31  }
32 }
33 
35 {
36  B2DEBUG(100, "In CDCMCMap::clear()");
37  m_simHitsByHit.clear();
38 
39  m_hitsByMCParticle.clear();
40  m_simHitsByMCParticle.clear();
41 
44 }
45 
47 {
48  B2DEBUG(100, "In CDCMCMap::fill()");
49  clear();
50 
54 
57 
58  B2DEBUG(100, "m_simHitsByHit.size(): " << m_simHitsByHit.size());
59 
60  B2DEBUG(100, "m_hitsByMCParticle.size(): " << m_hitsByMCParticle.size());
61  B2DEBUG(100, "m_simHitsByMCParticle.size(): " << m_simHitsByMCParticle.size());
62 
63  B2DEBUG(100, "m_reassignedSecondaryHits.size(): " << m_reassignedSecondaryHits.size());
64  B2DEBUG(100, "m_reassignedSecondarySimHits.size(): " << m_reassignedSecondarySimHits.size());
65 }
66 
68 {
69  StoreArray<CDCSimHit> simHits;
70  StoreArray<CDCHit> hits;
71 
72  // Pickup an iterator for hinted insertion
73  auto itInsertHint = m_simHitsByHit.end();
74 
75  for (const CDCSimHit& simHit : simHits) {
76  const CDCSimHit* ptrSimHit = &simHit;
77  RelationVector<CDCHit> relatedHits = simHit.getRelationsTo<CDCHit>();
78 
79  int nRelatedHits = relatedHits.size();
80  if (nRelatedHits > 1) {
81  B2WARNING("CDCSimHit as more than one related CDCHit - reorganize the mapping");
82  }
83 
84  for (const CDCHit& hit : relatedHits) {
85  const CDCHit* ptrHit = &hit;
86 
87  if (m_simHitsByHit.count(ptrHit) != 0) {
88  B2WARNING("CDCHit as more than one related CDCSimHit - reorganize the mapping");
89  }
90 
91  itInsertHint = m_simHitsByHit.insert(itInsertHint, {ptrHit, ptrSimHit});
92  }
93  }
94 
95  // Check if every hit has a corresponding simulated hit
96  for (const CDCHit& hit : hits) {
97  const CDCHit* ptrHit = &hit;
98 
99  if (m_simHitsByHit.count(ptrHit) == 0) {
100  B2WARNING("CDCHit has no related CDCSimHit in CDCMCMap::fill()");
101  }
102  }
103 }
104 
106 {
107  StoreArray<MCParticle> mcParticles;
108  StoreArray<CDCHit> hits;
109 
110  std::map<const MCParticle*, std::vector<const CDCSimHit*>> primarySimHitsByMCParticle;
111 
112  for (const CDCHit& hit : hits) {
113  const CDCHit* ptrHit = &hit;
114  RelationVector<MCParticle> relatedMCParticles = hit.getRelationsFrom<MCParticle>();
115 
116  const int nRelatedMCParticles = relatedMCParticles.size();
117 
118  if (nRelatedMCParticles == 0 and not isBackground(ptrHit)) {
119  B2WARNING("CDCHit has no related MCParticle but CDCHit indicates that it is no "
120  "background in CDCMCMap::fill()");
121  }
122 
123  if (nRelatedMCParticles > 1) {
124  B2WARNING("CDCHit as more than one related MCParticle - reorganize the mapping");
125  }
126 
127  // Pickup an iterator for hinted insertion
128  auto itInsertHint = m_hitsByMCParticle.end();
129  for (int iRelatedMCParticle = 0; iRelatedMCParticle < nRelatedMCParticles; ++iRelatedMCParticle) {
130  const MCParticle* ptrMCParticle = relatedMCParticles[iRelatedMCParticle];
131  double weight = relatedMCParticles.weight(iRelatedMCParticle);
132 
133  if (indicatesReassignedSecondary(weight)) {
134  m_reassignedSecondaryHits.insert(ptrHit);
135  } else {
136  const CDCSimHit* ptrSimHit = ptrHit->getRelatedFrom<CDCSimHit>();
137  if (ptrMCParticle->isPrimaryParticle() and ptrSimHit) {
138  primarySimHitsByMCParticle[ptrMCParticle].push_back(ptrSimHit);
139  }
140  }
141 
142  itInsertHint = m_hitsByMCParticle.insert(itInsertHint, {ptrMCParticle, ptrHit});
143  }
144  }
145 
146  // Check time ordering of primary hits
147  int nSortedIncorretly = 0;
148  auto lessFlightTime = [](const CDCSimHit * lhs, const CDCSimHit * rhs) {
149  return lhs->getFlightTime() < rhs->getFlightTime();
150  };
151 
152  auto lessArrayIndex = [](const CDCSimHit * lhs, const CDCSimHit * rhs) -> bool {
153  return lhs->getArrayIndex() < rhs->getArrayIndex();
154  };
155 
156  for (std::pair<const MCParticle* const, std::vector<const CDCSimHit*>>&
157  primarySimHitsForMCParticle : primarySimHitsByMCParticle) {
158 
159  const MCParticle* ptrMCParticle = primarySimHitsForMCParticle.first;
160  std::vector<const CDCSimHit*>& simHits = primarySimHitsForMCParticle.second;
161  std::sort(simHits.begin(), simHits.end(), lessArrayIndex);
162  auto itSorted = std::is_sorted_until(simHits.begin(), simHits.end(), lessFlightTime);
163  if (itSorted != simHits.end()) {
164  ++nSortedIncorretly;
165  B2DEBUG(100,
166  "CDCSimHits for MCParticle " << ptrMCParticle->getArrayIndex()
167  << " only sorted correctly up to hit number "
168  << std::distance(simHits.begin(), itSorted));
169  --itSorted;
170  B2DEBUG(100,
171  "Between wire " << (*itSorted)->getWireID() << " " << (*itSorted)->getFlightTime()
172  << "ns "
173  << (*itSorted)->getArrayIndex());
174  ++itSorted;
175  B2DEBUG(100,
176  "and wire " << (*itSorted)->getWireID() << " " << (*itSorted)->getFlightTime()
177  << "ns "
178  << (*itSorted)->getArrayIndex());
179  }
180  }
181  if (nSortedIncorretly) {
182  B2WARNING("(BII-2136) CDCSimHits for "
183  << nSortedIncorretly
184  << " primary mc particles are not sorted correctly by their time of flight");
185  }
186 }
187 
189 {
190  StoreArray<MCParticle> mcParticles;
191  StoreArray<CDCSimHit> simHits;
192 
193  for (const CDCSimHit& simHit : simHits) {
194  const CDCSimHit* ptrSimHit = &simHit;
195  RelationVector<MCParticle> relatedMCParticles = simHit.getRelationsFrom<MCParticle>();
196 
197  const int nRelatedMCParticles = relatedMCParticles.size();
198 
199  if (nRelatedMCParticles == 0 and not isBackground(ptrSimHit)) {
200  B2WARNING("CDCSimHit has no related MCParticle but CDCSimHit indicates that it is no "
201  "background in CDCMCMap::fill()");
202  }
203 
204  if (nRelatedMCParticles > 1) {
205  B2WARNING("CDCSimHit as more than one related MCParticle - reorganize the mapping");
206  }
207 
208  // Pickup an iterator for hinted insertion
209  auto itInsertHint = m_simHitsByMCParticle.end();
210  for (int iRelatedMCParticle = 0; iRelatedMCParticle < nRelatedMCParticles; ++iRelatedMCParticle) {
211  const MCParticle* ptrMCParticle = relatedMCParticles[iRelatedMCParticle];
212  double weight = relatedMCParticles.weight(iRelatedMCParticle);
213 
214  if (indicatesReassignedSecondary(weight)) {
215  m_reassignedSecondarySimHits.insert(ptrSimHit);
216  }
217 
218  itInsertHint = m_simHitsByMCParticle.insert(itInsertHint, {ptrMCParticle, ptrSimHit});
219  }
220  }
221 }
222 
224 {
225  StoreArray<CDCHit> hits;
226 
227  for (const CDCHit& hit : hits) {
228  const CDCHit* ptrHit = &hit;
229 
230  const CDCSimHit* ptrSimHit = getSimHit(ptrHit);
231  const MCParticle* ptrMCParticle = getMCParticle(ptrHit);
232 
233  const MCParticle* ptrMCParticleFromSimHit = getMCParticle(ptrSimHit);
234 
235  if (ptrMCParticle != ptrMCParticleFromSimHit) {
236  B2WARNING("MCParticle from CDCHit and MCParticle from related CDCSimHit mismatch in "
237  "CDCMCMap::validateRelations()");
238  }
239  }
240 }
241 
243 {
244  for (const CDCHit* ptrHit : m_reassignedSecondaryHits) {
245 
246  const CDCSimHit* ptrSimHit = getSimHit(ptrHit);
247  if (not isReassignedSecondary(ptrSimHit)) {
248  B2WARNING("CDCHit is reassigned secondary but related CDCSimHit is not.");
249  }
250  }
251 }
252 
253 MayBePtr<const CDCSimHit> CDCMCMap::getSimHit(const CDCHit* hit) const
254 {
255  return hit ? hit->getRelated<CDCSimHit>() : nullptr;
256 }
257 
258 MayBePtr<const CDCHit> CDCMCMap::getHit(const CDCSimHit* simHit) const
259 {
260  return simHit ? simHit->getRelated<CDCHit>() : nullptr;
261 }
262 
263 bool CDCMCMap::isBackground(const CDCSimHit* simHit) const
264 {
265  return simHit ? simHit->getBackgroundTag() != BackgroundMetaData::bg_none : false;
266 }
267 
268 bool CDCMCMap::isBackground(const CDCHit* hit) const
269 {
270  return isBackground(getSimHit(hit));
271 }
272 
273 MayBePtr<const MCParticle> CDCMCMap::getMCParticle(const CDCHit* hit) const
274 {
275  return hit ? hit->getRelated<MCParticle>() : nullptr;
276 }
277 
278 MayBePtr<const MCParticle> CDCMCMap::getMCParticle(const CDCSimHit* simHit) const
279 {
280  return simHit ? simHit->getRelated<MCParticle>() : nullptr;
281 }
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
Belle2::CDCSimHit::getFlightTime
double getFlightTime() const
The method to get flight time.
Definition: CDCSimHit.h:201
Belle2::TrackFindingCDC::CDCMCMap::m_simHitsByMCParticle
std::multimap< const MCParticle *, const CDCSimHit * > m_simHitsByMCParticle
Memory for a one to n relation from MCParticles to CDCSimHit.
Definition: CDCMCMap.h:169
Belle2::TrackFindingCDC::CDCMCMap::fill
void fill()
Fill the Monte Carlo information retrieved from the DataStore into the local multimaps.
Definition: CDCMCMap.cc:46
Belle2::TrackFindingCDC::CDCMCMap::getHit
MayBePtr< const CDCHit > getHit(const CDCSimHit *simHit) const
Seeks the CDCHit related to the CDCSimHit - nullptr if no CDCHit is related.
Definition: CDCMCMap.cc:258
Belle2::CDCHit
Class containing the result of the unpacker in raw data and the result of the digitizer in simulation...
Definition: CDCHit.h:51
Belle2::TrackFindingCDC::CDCMCMap::m_hitsByMCParticle
std::multimap< const MCParticle *, const CDCHit * > m_hitsByMCParticle
Memory for a one to n relation from MCParticles to CDCHit.
Definition: CDCMCMap.h:166
Belle2::TrackFindingCDC::CDCMCMap::m_simHitsByHit
std::multimap< const CDCHit *, const CDCSimHit * > m_simHitsByHit
Memory for a one to one relation from CDCHit to CDCSimHits.
Definition: CDCMCMap.h:163
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:280
Belle2::TrackFindingCDC::CDCMCMap::getMCParticle
MayBePtr< const MCParticle > getMCParticle(const CDCHit *hit) const
Seeks the MCParticle related to the CDCHit.
Definition: CDCMCMap.cc:273
Belle2::CDCSimHit
Example Detector.
Definition: CDCSimHit.h:33
Belle2::MCParticle::getArrayIndex
int getArrayIndex() const
Get 0-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:255
Belle2::MCParticle::isPrimaryParticle
bool isPrimaryParticle() const
Check if particle is a primary particle which was created by the generator (and not,...
Definition: MCParticle.h:588
Belle2::TrackFindingCDC::CDCMCMap::indicatesReassignedSecondary
static bool indicatesReassignedSecondary(double weight)
Indicate if the given weight suggests that the corresponding hit to MCParticle relation has been redi...
Definition: CDCMCMap.h:62
Belle2::TrackFindingCDC::CDCMCMap::m_reassignedSecondarySimHits
std::set< const CDCSimHit * > m_reassignedSecondarySimHits
The set of reassigned secondary CDCSimHits.
Definition: CDCMCMap.h:175
Belle2::TrackFindingCDC::CDCMCMap::fillMCParticleBySimHitMap
void fillMCParticleBySimHitMap()
Retrieve the relations array from MCParticle to CDCSimHits and fill it in to the local map which does...
Definition: CDCMCMap.cc:188
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::CDCMCMap::validateReassignedSecondaries
void validateReassignedSecondaries() const
Checks if each CDCHit is marked as reassigned secondary is related to a reassigned secondary CDCSimHi...
Definition: CDCMCMap.cc:242
Belle2::TrackFindingCDC::CDCMCMap::validateRelations
void validateRelations() const
Checks if the relations CDCHit -> MCParticle and CDCHit -> CDCSimHit -> MCParticle commute.
Definition: CDCMCMap.cc:223
Belle2::TrackFindingCDC::CDCMCMap::m_reassignedSecondaryHits
std::set< const CDCHit * > m_reassignedSecondaryHits
The set of reassigned secondary CDCHits.
Definition: CDCMCMap.h:172
Belle2::RelationsInterface::getArrayIndex
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Definition: RelationsObject.h:387
Belle2::BackgroundMetaData::bg_none
@ bg_none
No background.
Definition: BackgroundMetaData.h:41
Belle2::TrackFindingCDC::CDCMCMap::fillMCParticleByHitMap
void fillMCParticleByHitMap()
Retrieve the relations array from MCParticle to CDCHits and fill it in to the local map which does th...
Definition: CDCMCMap.cc:105
Belle2::TrackFindingCDC::CDCMCMap::getSimHit
MayBePtr< const CDCSimHit > getSimHit(const CDCHit *hit) const
Seeks the CDCSimHit related to the CDCHit.
Definition: CDCMCMap.cc:253
Belle2::TrackFindingCDC::CDCMCMap::clear
void clear()
Clear all information from the former event.
Definition: CDCMCMap.cc:34
Belle2::RelationVector::weight
float weight(int index) const
Get weight with index.
Definition: RelationVector.h:120
Belle2::SimHitBase::getBackgroundTag
virtual unsigned short getBackgroundTag() const
Get background tag.
Definition: SimHitBase.h:56
Belle2::TrackFindingCDC::CDCMCMap::isReassignedSecondary
bool isReassignedSecondary(const CDCSimHit *ptrSimHit) const
Indicates if the CDCSimHit has been reassigned to a primary MCParticle.
Definition: CDCMCMap.h:120
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::TrackFindingCDC::CDCMCMap::isBackground
bool isBackground(const CDCSimHit *simHit) const
Indicates if the CDCSimHit is considered background.
Definition: CDCMCMap.cc:263
Belle2::TrackFindingCDC::CDCMCMap::fillSimHitByHitMap
void fillSimHitByHitMap()
Retrieve the relations array from CDCSimHits to CDCHits and fill it in to the local map which does th...
Definition: CDCMCMap.cc:67
Belle2::RelationsInterface::getRelatedFrom
FROM * getRelatedFrom(const std::string &name="", const std::string &namedRelation="") const
Get the object from which this object has a relation.
Definition: RelationsObject.h:265