11 #include <trg/klm/modules/klmtrigger/KLMTriggerModule.h>
14 #include <framework/datastore/StoreArray.h>
15 #include <framework/datastore/StoreObjPtr.h>
18 #include <framework/dataobjects/EventMetaData.h>
21 #include <klm/dataobjects/KLMDigit.h>
23 #include <trg/klm/dataobjects/KLMTriggerHit.h>
24 #include <trg/klm/dataobjects/KLMTriggerTrack.h>
26 #include <unordered_map>
41 KLMTriggerModule::KLMTriggerModule() :
Module()
46 "Maximum chi squared for a track",
49 "Maximum impact parameter for a track",
52 "Minimum number of fired layers for a track",
63 klmTriggerHits.registerInDataStore();
67 klmTriggerTracks.registerInDataStore();
74 B2DEBUG(100,
"KLMTrigger: Experiment " << evtMetaData->getExperiment() <<
", run " << evtMetaData->getRun());
102 for (
int i = 0; i < nEntries; ++i) {
103 const KLMDigit* bklmDigit_i = klmDigits[i];
106 for (
int j = i + 1; j < nEntries; ++j) {
107 const KLMDigit* bklmDigit_j = klmDigits[j];
116 int sector = bklmDigit_i->
getSector() - 1;
117 int layer = bklmDigit_i->
getLayer() - 1;
121 phiStrip = bklmDigit_i->
getStrip() - 1;
122 zStrip = bklmDigit_j->
getStrip() - 1;
124 zStrip = bklmDigit_i->
getStrip() - 1;
125 phiStrip = bklmDigit_j->
getStrip() - 1;
128 int xInt = 0, yInt = 0, zInt = 0;
129 double x = 0.0, y = 0.0, z = 0.0;
132 x = (double)(xInt >> 3);
133 y = (double)(yInt >> 3);
134 z = (double)(zInt >> 3);
143 hit->addRelationTo(klmDigits[i]);
144 hit->addRelationTo(klmDigits[j]);
162 std::unordered_map<int, KLMTriggerTrack*> trackMap;
166 for (
int i = 0; i < nEntries; ++i) {
169 int section = hit->getSection();
170 int sector = hit->getSector();
174 if (trackMap.find(sectorID) == trackMap.end())
175 trackMap[sectorID] = klmTriggerTracks.
appendNew(section, sector);
177 trackMap[sectorID]->addRelationTo(klmTriggerHits[i]);
185 if (!klmTriggerTracks.
isValid())
189 for (
int i = 0; i < nEntries; ++i) {
193 int nHits = hits.size();
195 int sumX = 0, sumY = 0, sumZ = 0, sumXX = 0, sumXY = 0, sumXZ = 0, sumYY = 0, sumZZ = 0;
200 std::unordered_map<int, bool> layersMap;
202 for (
int j = 0; j < nHits; ++j) {
205 const int xInt = hit->getXInt();
206 const int yInt = hit->getYInt();
207 const int zInt = hit->getZInt();
208 const int layer = hit->getLayer();
214 sumXX += xInt * xInt;
215 sumXY += xInt * yInt;
216 sumXZ += xInt * zInt;
217 sumYY += yInt * yInt;
218 sumZZ += zInt * zInt;
221 if (layer < firstLayer)
223 if (layer > lastLayer)
226 layersMap[layer] =
true;
259 int denom = sumXX * nHits - sumX * sumX;
261 double denomInversed = 0.0;
262 double slopeXY = 0.0;
263 double interceptXY = 0.0;
264 double chisqXY = 0.0;
266 double slopeXZ = 0.0;
267 double interceptXZ = 0.0;
268 double chisqXZ = 0.0;
274 denomInversed = 1.0 / denom;
276 slopeXY = denomInversed * (sumXY * nHits - sumX * sumY);
277 interceptXY = denomInversed * (sumXX * sumY - sumX * sumXY);
279 slopeXZ = denomInversed * (sumXZ * nHits - sumX * sumZ);
280 interceptXZ = denomInversed * (sumXX * sumZ - sumX * sumXZ);
293 ipXY = interceptXY * interceptXY * (1.0 - slopeXY * slopeXY);
294 chisqXY = slopeXY * slopeXY * sumXX + interceptXY * interceptXY * nHits + sumYY + 2.0 * slopeXY * interceptXY * sumX - 2.0 *
295 slopeXY * sumXY - 2.0 * interceptXY * sumY;
297 ipXZ = interceptXZ * interceptXZ * (1.0 - slopeXZ * slopeXZ);
298 chisqXZ = slopeXZ * slopeXZ * sumXX + interceptXZ * interceptXZ * nHits + sumZZ + 2.0 * slopeXZ * interceptXZ * sumX - 2.0 *
299 slopeXZ * sumXZ - 2.0 * interceptXZ * sumZ;
303 interceptXY *= (32 / 8);
304 chisqXY *= (1024 / 64);
307 interceptXZ *= (32 / 8);
308 chisqXZ *= (1024 / 64);
316 nLayers = layersMap.size();
318 track->setSlopeXY(slopeXY);
319 track->setInterceptXY(interceptXY);
320 track->setChisqXY(chisqXY);
321 track->setImpactParameterXY(ipXY);
323 track->setSlopeXZ(slopeXZ);
324 track->setInterceptXZ(interceptXZ);
325 track->setChisqXZ(chisqXZ);
326 track->setImpactParameterXZ(ipXZ);
328 track->setFirstLayer(firstLayer);
329 track->setLastLayer(lastLayer);
330 track->setNLayers(nLayers);
338 track->setTrigger(
true);
409 const int c_LayerXCoord[
c_TotalLayers] = {1628, 1700, 1773, 1846, 1919, 1992, 2064, 2137, 2210, 2283, 2356, 2428, 2501, 2574, 2647};
410 const int c_LayerY0[
c_TotalLayers] = { -2403, -2566, -2744, -2862, -2979, -3097, -3234, -3354, -3474, -3587, -3708, -3828, -3948, -4061, -4181};
411 const int c_PhiWidth[
c_TotalLayers] = {128, 128, 157, 164, 170, 177, 138, 143, 148, 153, 158, 163, 168, 173, 178};
413 const int c_Z01 = 18;
414 const int c_Z02 = 25;
415 const int c_ZWidth1 = 32;
416 const int c_ZWidth2 = 36;
417 const int c_ZOffset = 376;
419 bool flipped =
false;
426 if (layer == 0 && section == 0 && sector == 0)
428 else if (layer == 0 && section == 0 && sector == 1)
430 else if (layer == 0 && section == 0 && sector == 2)
432 else if (layer == 0 && section == 0 && sector == 3)
434 else if (layer == 0 && section == 0 && sector == 4)
436 else if (layer == 0 && section == 0 && sector == 5)
438 else if (layer == 0 && section == 0 && sector == 6)
440 else if (layer == 0 && section == 0 && sector == 7)
442 else if (layer == 0 && section == 1 && sector == 0)
444 else if (layer == 0 && section == 1 && sector == 1)
446 else if (layer == 0 && section == 1 && sector == 2)
448 else if (layer == 0 && section == 1 && sector == 3)
450 else if (layer == 0 && section == 1 && sector == 4)
452 else if (layer == 0 && section == 1 && sector == 5)
454 else if (layer == 0 && section == 1 && sector == 6)
456 else if (layer == 0 && section == 1 && sector == 7)
458 else if (layer == 1 && section == 0 && sector == 0)
460 else if (layer == 1 && section == 0 && sector == 1)
462 else if (layer == 1 && section == 0 && sector == 2)
464 else if (layer == 1 && section == 0 && sector == 3)
466 else if (layer == 1 && section == 0 && sector == 4)
468 else if (layer == 1 && section == 0 && sector == 5)
470 else if (layer == 1 && section == 0 && sector == 6)
472 else if (layer == 1 && section == 0 && sector == 7)
474 else if (layer == 1 && section == 1 && sector == 0)
476 else if (layer == 1 && section == 1 && sector == 1)
478 else if (layer == 1 && section == 1 && sector == 2)
480 else if (layer == 1 && section == 1 && sector == 3)
482 else if (layer == 1 && section == 1 && sector == 4)
484 else if (layer == 1 && section == 1 && sector == 5)
486 else if (layer == 1 && section == 1 && sector == 6)
488 else if (layer == 1 && section == 1 && sector == 7)
490 else if (layer > 2 && section == 0)
505 if (layer == 2 || layer == 4 || layer == 8 ||
506 layer == 10 || layer == 11 || layer == 13) {
509 else if (phiStrip > 40)
511 else if (phiStrip > 35)
513 else if (phiStrip > 29)
515 else if (phiStrip > 24)
517 else if (phiStrip > 18)
519 else if (phiStrip > 13)
521 else if (phiStrip > 8)
523 else if (phiStrip > 2)
527 }
else if (layer == 3 || layer == 6 or
528 layer == 7 || layer == 9) {
531 else if (phiStrip > 45)
533 else if (phiStrip > 42)
535 else if (phiStrip > 39)
537 else if (phiStrip > 36)
539 else if (phiStrip > 34)
541 else if (phiStrip > 31)
543 else if (phiStrip > 28)
545 else if (phiStrip > 25)
547 else if (phiStrip > 23)
549 else if (phiStrip > 20)
551 else if (phiStrip > 17)
553 else if (phiStrip > 15)
555 else if (phiStrip > 12)
557 else if (phiStrip > 9)
559 else if (phiStrip > 6)
561 else if (phiStrip > 4)
563 else if (phiStrip > 1)
576 else if (zStrip > 38)
578 else if (zStrip > 32)
580 else if (zStrip > 26)
582 else if (zStrip > 20)
584 else if (zStrip > 13)
595 y = c_LayerY0[layer] + dy + phiStrip * c_PhiWidth[layer];
599 z = z0 + dz + zStrip * zWidth;
605 x = c_LayerXCoord[layer];