Belle II Software development
KLMStripEfficiencyCollectorModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9/* Own header. */
10#include <klm/modules/KLMStripEfficiencyCollector/KLMStripEfficiencyCollectorModule.h>
11#include <klm/dataobjects/KLMChannelIndex.h>
12
13/* CLHEP headers. */
14#include <CLHEP/Vector/ThreeVector.h>
15
16using namespace Belle2;
17
18REG_MODULE(KLMStripEfficiencyCollector);
19
22 m_ElementNumbers(&(KLMElementNumbers::Instance())),
24 m_GeometryBKLM(nullptr),
26 m_MatchingFile(nullptr),
27 m_MatchingTree(nullptr),
28 m_MatchingHitData({0, 0, 0, 0, 0, 0, 0., nullptr, nullptr}),
29 m_MatchedStrip(0)
30{
31 setDescription("Module for KLM strip efficiency data collection.");
32 addParam("MuonListName", m_MuonListName, "Muon list name.",
33 std::string("mu+:all"));
34 addParam("AllowedDistance1D", m_AllowedDistance1D,
35 "Maximal distance in the units of strip number from ExtHit to "
36 "matching KLMDigit (not for multi-strip hits).", double(8));
37 addParam("MinimalMatchingDigits", m_MinimalMatchingDigits,
38 "Minimal number of matching digits.", 0);
39 addParam("MinimalMatchingDigitsOuterLayers",
40 m_MinimalMatchingDigitsOuterLayers,
41 "Minimal number of matching digits in outer layers.", 0);
42 addParam("MinimalMomentumNoOuterLayers", m_MinimalMomentumNoOuterLayers,
43 "Minimal momentum in case there are no hits in outer layers.", 0.0);
44 addParam("RemoveUnusedMuons", m_RemoveUnusedMuons,
45 "Whether to remove unused muons.", false);
46 addParam("IgnoreBackwardPropagation", m_IgnoreBackwardPropagation,
47 "Whether to ignore ExtHits with backward propagation.", false);
48 addParam("Debug", m_Debug, "Debug mode.", false);
49 addParam("DebugFileName", m_MatchingFileName, "Debug file name.", std::string("matching.root"));
50 setPropertyFlags(c_ParallelProcessingCertified);
51}
52
56
58{
59 m_Digits.isRequired();
60 m_tracks.isRequired();
61 m_extHits.isRequired();
62 m_MuonList.isRequired(m_MuonListName);
63 int nPlanes = m_PlaneArrayIndex->getNElements();
64 TH1F* matchedDigitsInPlane = new TH1F(
65 "matchedDigitsInPlane", "Number of matching KLMDigits",
66 nPlanes, -0.5, double(nPlanes) - 0.5);
67 TH1F* allExtHitsInPlane = new TH1F(
68 "allExtHitsInPlane", "Number of ExtHits",
69 nPlanes, -0.5, double(nPlanes) - 0.5);
70 registerObject<TH1F>("matchedDigitsInPlane", matchedDigitsInPlane);
71 registerObject<TH1F>("allExtHitsInPlane", allExtHitsInPlane);
73 if (m_Debug) {
74 m_MatchingFile = new TFile(m_MatchingFileName.c_str(), "recreate");
75 m_MatchingTree = new TTree("t_matching", "");
76 m_MatchingTree->Branch("subdetector", &m_MatchingHitData.subdetector,
77 "subdetector/I");
78 m_MatchingTree->Branch("section", &m_MatchingHitData.section, "section/I");
79 m_MatchingTree->Branch("sector", &m_MatchingHitData.sector, "sector/I");
80 m_MatchingTree->Branch("layer", &m_MatchingHitData.layer, "layer/I");
81 m_MatchingTree->Branch("plane", &m_MatchingHitData.plane, "plane/I");
82 m_MatchingTree->Branch("strip", &m_MatchingHitData.strip, "strip/I");
83 m_MatchingTree->Branch("matchedStrip", &m_MatchedStrip, "matchedStrip/I");
84 }
85}
86
88{
89 if (m_Debug) {
90 m_MatchingFile->cd();
91 m_MatchingTree->Write();
92 delete m_MatchingTree;
93 delete m_MatchingFile;
94 }
95}
96
98{
99 if (!m_ChannelStatus.isValid())
100 B2FATAL("KLM channel status data are not available.");
101 int minimalActivePlanes = -1;
103 for (KLMChannelIndex& klmSector : klmSectors) {
104 KLMChannelIndex klmNextSector(klmSector);
105 ++klmNextSector;
106 int activePlanes = 0;
107 KLMChannelIndex klmPlane(klmSector);
109 for (; klmPlane != klmNextSector; ++klmPlane) {
110 KLMChannelIndex klmNextPlane(klmPlane);
111 ++klmNextPlane;
112 bool isActive = false;
113 KLMChannelIndex klmChannel(klmPlane);
115 for (; klmChannel != klmNextPlane; ++ klmChannel) {
116 KLMChannelNumber channel = klmChannel.getKLMChannelNumber();
118 m_ChannelStatus->getChannelStatus(channel);
119 if (status == KLMChannelStatus::c_Unknown)
120 B2FATAL("Incomplete KLM channel status data.");
121 if (status == KLMChannelStatus::c_Normal) {
122 isActive = true;
123 break;
124 }
125 }
126 if (isActive)
127 ++activePlanes;
128 }
129 /*
130 * If a sector is completely off, it is not necessary to reduce
131 * the minimal number of digits, because the efficiencies cannot be
132 * measured for the entire sector anyway.
133 */
134 if (activePlanes == 0)
135 continue;
136 if (minimalActivePlanes < 0)
137 minimalActivePlanes = activePlanes;
138 else if (minimalActivePlanes < activePlanes)
139 minimalActivePlanes = activePlanes;
140 }
141 if ((minimalActivePlanes >= 0) &&
142 (minimalActivePlanes < m_MinimalMatchingDigits)) {
143 B2WARNING("The minimal number of active planes (" << minimalActivePlanes <<
144 ") is less than the minimal number of matching digits (" <<
145 m_MinimalMatchingDigits << "). The minimal number of "
146 "matching digits is reduced to make the calibration possible.");
147 m_MinimalMatchingDigits = minimalActivePlanes;
148 }
149}
150
154
156{
157 std::vector<unsigned int> toRemove;
158 unsigned int nMuons = m_MuonList->getListSize();
159 /*
160 * The getter functions create the object for particular experiment and run.
161 * It is not guaranteed that collectDataTrack() is called if there are
162 * no tracks (e.g. for too short or bad runs). Thus, they are called here
163 * to guarantee the creation of histograms.
164 */
165 TH1F* matchedDigitsInPlane = getObjectPtr<TH1F>("matchedDigitsInPlane");
166 TH1F* allExtHitsInPlane = getObjectPtr<TH1F>("allExtHitsInPlane");
167 for (unsigned int i = 0; i < nMuons; ++i) {
168 const Particle* muon = m_MuonList->getParticle(i);
169 bool trackUsed = collectDataTrack(muon, matchedDigitsInPlane,
170 allExtHitsInPlane);
171 if (m_RemoveUnusedMuons && !trackUsed)
172 toRemove.push_back(muon->getArrayIndex());
173 }
175 m_MuonList->removeParticles(toRemove);
176}
177
179 std::map<KLMPlaneNumber, struct HitData>& hitMap,
180 KLMPlaneNumber planeGlobal, struct HitData* hitData)
181{
182 std::map<KLMPlaneNumber, struct HitData>::iterator it;
183 it = hitMap.find(planeGlobal);
184 /*
185 * There may be more than one such hit e.g. if track crosses the edge
186 * of the strips or WLS fiber groove. Select only one hit per plane.
187 */
188 if (it == hitMap.end()) {
189 hitMap.insert(std::pair<KLMPlaneNumber, struct HitData>(
190 planeGlobal, *hitData));
191 }
192}
193
195 struct HitData* hitData)
196{
197 for (const KLMDigit& digit : m_Digits) {
198 if (!(digit.getSubdetector() == hitData->subdetector &&
199 digit.getSection() == hitData->section &&
200 digit.getLayer() == hitData->layer &&
201 digit.getSector() == hitData->sector &&
202 digit.getPlane() == hitData->plane))
203 continue;
204
205 double stripPosition = digit.getStrip();
206 double allowedDistance1D = m_AllowedDistance1D;
207
208 if (digit.isMultiStrip()) {
209 stripPosition = 0.5 * (digit.getLastStrip() + digit.getStrip());
210 allowedDistance1D *= (digit.getLastStrip() - digit.getStrip() + 1);
211 }
212 if (fabs(stripPosition - hitData->strip) < allowedDistance1D) {
213 hitData->digit = &digit;
214 return;
215 }
216 }
217}
218
220 const Particle* muon, TH1F* matchedDigitsInPlane, TH1F* allExtHitsInPlane)
221{
222 const int nExtrapolationLayers =
224 const Track* track = muon->getTrack();
225 RelationVector<ExtHit> extHits = track->getRelationsTo<ExtHit>();
226 std::map<KLMPlaneNumber, struct HitData> selectedHits;
227 std::map<KLMPlaneNumber, struct HitData>::iterator it;
228 KLMChannelNumber channel;
230 struct HitData hitData, hitDataPrevious;
231 ROOT::Math::XYZVector extHitPosition;
232 CLHEP::Hep3Vector extHitPositionCLHEP, localPosition;
233 int layer;
234 int extHitLayer[nExtrapolationLayers] = {0};
235 int digitLayer[nExtrapolationLayers] = {0};
236 // initialize hitDataPrevious components
237 hitDataPrevious.subdetector = -1;
238 hitDataPrevious.section = -1;
239 hitDataPrevious.sector = -1;
240 hitDataPrevious.layer = -1;
241 for (const ExtHit& hit : extHits) {
242 /*
243 * Choose hits that exit the sensitive volume.
244 * It is not possible to use entry hits because of a bug in Geant4:
245 * the step does not always have a correct status (fGeomBoundary),
246 * and, consequently, ExtHits are not created.
247 */
248 if (hit.getStatus() != EXT_EXIT)
249 continue;
250 /*
251 * Ignore ExtHits with backward propagation. This affects cosmic events
252 * only. The removal of hits with backward propagation is normally
253 * not needed, however, it is added because of backward error propagation
254 * bug in Geant4 10.6.
255 */
257 if (hit.isBackwardPropagated())
258 continue;
259 }
260 KLMPlaneNumber planeGlobal = 0;
261 hitData.hit = &hit;
262 hitData.digit = nullptr;
263 if (hit.getDetectorID() == Const::EDetector::EKLM) {
264 int stripGlobal = hit.getCopyID();
266 m_eklmElementNumbers->stripNumberToElementNumbers(
267 stripGlobal, &hitData.section, &hitData.layer, &hitData.sector,
268 &hitData.plane, &hitData.strip);
269 channel = m_ElementNumbers->channelNumberEKLM(
270 hitData.section, hitData.sector, hitData.layer,
271 hitData.plane, hitData.strip);
272 status = m_ChannelStatus->getChannelStatus(channel);
273 if (status == KLMChannelStatus::c_Unknown)
274 B2FATAL("Incomplete KLM channel status data.");
275 if (status == KLMChannelStatus::c_Normal) {
276 layer = m_ElementNumbers->getExtrapolationLayer(
277 hitData.subdetector, hitData.layer);
278 extHitLayer[layer - 1]++;
279 planeGlobal = m_ElementNumbers->planeNumberEKLM(
280 hitData.section, hitData.sector, hitData.layer,
281 hitData.plane);
282 addHit(selectedHits, planeGlobal, &hitData);
283 }
284 } else if (hit.getDetectorID() == Const::EDetector::BKLM) {
285 int moduleNumber = hit.getCopyID();
288 moduleNumber, &hitData.section, &hitData.sector, &hitData.layer);
290 /*
291 * For scintillators, the plane and strip numbers are recorded
292 * in the copy number.
293 */
295 moduleNumber, &hitData.section, &hitData.sector, &hitData.layer,
296 &hitData.plane, &hitData.strip);
297 channel = m_ElementNumbers->channelNumberBKLM(
298 hitData.section, hitData.sector, hitData.layer,
299 hitData.plane, hitData.strip);
300 status = m_ChannelStatus->getChannelStatus(channel);
301 if (status == KLMChannelStatus::c_Unknown)
302 B2FATAL("Incomplete KLM channel status data.");
303 if (status == KLMChannelStatus::c_Normal) {
304 layer = m_ElementNumbers->getExtrapolationLayer(
305 hitData.subdetector, hitData.layer);
306 extHitLayer[layer - 1]++;
307 planeGlobal = m_ElementNumbers->planeNumberBKLM(
308 hitData.section, hitData.sector, hitData.layer,
309 hitData.plane);
310 addHit(selectedHits, planeGlobal, &hitData);
311 }
312 } else {
313 /* For RPCs, the sensitive volume corresponds to both readout planes. */
314 extHitPosition = hit.getPosition();
315 extHitPositionCLHEP.setX(extHitPosition.X());
316 extHitPositionCLHEP.setY(extHitPosition.Y());
317 extHitPositionCLHEP.setZ(extHitPosition.Z());
318 const bklm::Module* module =
319 m_GeometryBKLM->findModule(hitData.section, hitData.sector,
320 hitData.layer);
321 localPosition = module->globalToLocal(extHitPositionCLHEP);
323 hitData.strip = module->getZStrip(localPosition);
324 /*
325 * FIXME
326 * There are 2 hits per module in RPCs, but the plane information is
327 * not available in ExtHit. For now, 2 entries are created (one for
328 * each plane) for the first hit, and the second one is removed.
329 */
330 if ((hitData.subdetector == hitDataPrevious.subdetector) &&
331 (hitData.section == hitDataPrevious.section) &&
332 (hitData.sector == hitDataPrevious.sector) &&
333 (hitData.layer == hitDataPrevious.layer))
334 continue;
335 std::memcpy(&hitDataPrevious, &hitData, sizeof(struct HitData));
336 /* The returned strip may be out of the valid range. */
338 hitData.section, hitData.sector, hitData.layer, hitData.plane,
339 hitData.strip, false)) {
340 channel = m_ElementNumbers->channelNumberBKLM(
341 hitData.section, hitData.sector, hitData.layer,
342 hitData.plane, hitData.strip);
343 status = m_ChannelStatus->getChannelStatus(channel);
344 if (status == KLMChannelStatus::c_Unknown)
345 B2FATAL("Incomplete KLM channel status data.");
346 if (status == KLMChannelStatus::c_Normal) {
347 layer = m_ElementNumbers->getExtrapolationLayer(
348 hitData.subdetector, hitData.layer);
349 extHitLayer[layer - 1]++;
350 hitData.localPosition = localPosition.z();
351 planeGlobal = m_ElementNumbers->planeNumberBKLM(
352 hitData.section, hitData.sector, hitData.layer,
353 hitData.plane);
354 addHit(selectedHits, planeGlobal, &hitData);
355 }
356 }
358 hitData.strip = module->getPhiStrip(localPosition);
359 /* The returned strip may be out of the valid range. */
361 hitData.section, hitData.sector, hitData.layer, hitData.plane,
362 hitData.strip, false)) {
363 channel = m_ElementNumbers->channelNumberBKLM(
364 hitData.section, hitData.sector, hitData.layer,
365 hitData.plane, hitData.strip);
366 status = m_ChannelStatus->getChannelStatus(channel);
367 if (status == KLMChannelStatus::c_Unknown)
368 B2FATAL("Incomplete KLM channel status data.");
369 if (status == KLMChannelStatus::c_Normal) {
370 layer = m_ElementNumbers->getExtrapolationLayer(
371 hitData.subdetector, hitData.layer);
372 extHitLayer[layer - 1]++;
373 hitData.localPosition = localPosition.y();
374 planeGlobal = m_ElementNumbers->planeNumberBKLM(
375 hitData.section, hitData.sector, hitData.layer,
376 hitData.plane);
377 addHit(selectedHits, planeGlobal, &hitData);
378 }
379 }
380 }
381 } else
382 continue;
383 }
384 /* Find matching digits. */
385 int nDigits = 0;
386 for (it = selectedHits.begin(); it != selectedHits.end(); ++it) {
387 findMatchingDigit(&(it->second));
388 if (it->second.digit != nullptr) {
389 nDigits++;
390 layer = m_ElementNumbers->getExtrapolationLayer(
391 it->second.subdetector, it->second.layer);
392 digitLayer[layer - 1]++;
393 }
394 }
395 if (nDigits < m_MinimalMatchingDigits)
396 return false;
397 /* Write efficiency histograms */
398 for (it = selectedHits.begin(); it != selectedHits.end(); ++it) {
399 int matchingDigits = 0;
400 int matchingDigitsOuterLayers = 0;
401 int extHitsOuterLayers = 0;
402 layer = m_ElementNumbers->getExtrapolationLayer(
403 it->second.subdetector, it->second.layer) - 1;
404 for (int i = 0; i < nExtrapolationLayers; ++i) {
405 if (i != layer)
406 matchingDigits += digitLayer[i];
407 if (i > layer) {
408 matchingDigitsOuterLayers += digitLayer[i];
409 extHitsOuterLayers += extHitLayer[i];
410 }
411 }
412 /* Check the number of matching digits in other layers. */
413 if (matchingDigits < m_MinimalMatchingDigits)
414 continue;
415 /*
416 * Check the number of matching digits in outer layers relatively to
417 * this hit.
418 */
419 if (matchingDigitsOuterLayers < m_MinimalMatchingDigitsOuterLayers) {
425 if (extHitsOuterLayers >= m_MinimalMatchingDigitsOuterLayers)
426 continue;
427 /*
428 * If the number of ExtHits is insufficient, then check the momentum.
429 * The muons with sufficiently large momentum have a very small
430 * probability to get absorbed in the detector.
431 */
433 continue;
434 }
435 allExtHitsInPlane->Fill(m_PlaneArrayIndex->getIndex(it->first));
436 if (it->second.digit) {
437 matchedDigitsInPlane->Fill(m_PlaneArrayIndex->getIndex(it->first));
438 if (m_Debug) {
439 std::memcpy(&m_MatchingHitData, &(it->second), sizeof(struct HitData));
440 m_MatchedStrip = it->second.digit->getStrip();
441 m_MatchingTree->Fill();
442 }
443 }
444 }
445 return true;
446}
static void channelNumberToElementNumbers(KLMChannelNumber channel, int *section, int *sector, int *layer, int *plane, int *strip)
Get element numbers by channel number.
static bool checkChannelNumber(int section, int sector, int layer, int plane, int strip, bool fatalError=true)
Check channel number.
@ c_FirstRPCLayer
First RPC layer.
static void moduleNumberToElementNumbers(KLMModuleNumber module, int *section, int *sector, int *layer)
Get element numbers by module number.
void registerObject(std::string name, T *obj)
Register object with a name, takes ownership, do not access the pointer beyond prepare()
CalibrationCollectorModule()
Constructor. Sets the default prefix for calibration dataobjects.
T * getObjectPtr(std::string name)
Calls the CalibObjManager to get the requested stored collector data.
Store one Ext hit as a ROOT object.
Definition ExtHit.h:31
KLM channel index.
KLMChannelNumber getKLMChannelNumber() const
Get KLM channel number.
void setIndexLevel(enum IndexLevel indexLevel)
Set index level.
@ c_Normal
Normally operating channel.
@ c_Unknown
Unknown status (no data).
KLM digit (class representing a digitized hit in RPCs or scintillators).
Definition KLMDigit.h:29
static constexpr int getMaximalExtrapolationLayer()
Get maximal extrapolation layer.
KLM plane array index.
const bklm::GeometryPar * m_GeometryBKLM
BKLM geometry.
void addHit(std::map< KLMPlaneNumber, struct HitData > &hitMap, KLMPlaneNumber planeGlobal, struct HitData *hitData)
Add hit to map.
DBObjPtr< KLMChannelStatus > m_ChannelStatus
Channel status.
void startRun() override
This method is called at the beginning of the run.
double m_MinimalMomentumNoOuterLayers
Minimal momentum in case there are no hits in outer layers.
const KLMElementNumbers * m_ElementNumbers
Element numbers.
const KLMPlaneArrayIndex * m_PlaneArrayIndex
Plane array index.
const EKLMElementNumbers * m_eklmElementNumbers
EKLM element numbers.
void findMatchingDigit(struct HitData *hitData)
Find matching digit.
void collect() override
This method is called for each event.
bool collectDataTrack(const Particle *muon, TH1F *matchedDigitsInPlane, TH1F *allExtHitsInPlane)
Collect the data for one muon.
void closeRun() override
This method is called at the end of run.
bool m_IgnoreBackwardPropagation
Whether to ignore ExtHits with backward propagation.
double m_AllowedDistance1D
Maximal distance in the units of strip number from ExtHit to matching KLMDigit.
int m_MinimalMatchingDigitsOuterLayers
Minimal number of matching digits in outer layers.
int m_MinimalMatchingDigits
Minimal number of matching digits.
Class to store reconstructed particles.
Definition Particle.h:76
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition Particle.cc:916
double getP() const
Returns momentum magnitude (same as getMomentumMagnitude but with shorter name)
Definition Particle.h:598
Class for type safe access to objects that are referred to in relations.
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Class that bundles various TrackFitResults.
Definition Track.h:25
static GeometryPar * instance(void)
Static method to get a reference to the singleton GeometryPar instance.
Define the geometry of a BKLM module Each sector [octant] contains Modules.
Definition Module.h:76
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
uint16_t KLMChannelNumber
Channel number.
uint16_t KLMPlaneNumber
Plane number.
Abstract base class for different kinds of events.