|
19 | 19 | #include "PWGUD/Core/DGCutparHolder.h" |
20 | 20 | #include "PWGUD/Core/UPCHelpers.h" |
21 | 21 |
|
| 22 | +#include "Common/Core/RecoDecay.h" |
22 | 23 | #include "Common/DataModel/EventSelection.h" |
23 | 24 | #include "Common/DataModel/TrackSelectionTables.h" |
24 | 25 |
|
|
29 | 30 |
|
30 | 31 | #include "TLorentzVector.h" |
31 | 32 |
|
| 33 | +#include <algorithm> |
32 | 34 | #include <bitset> |
33 | 35 | #include <vector> |
34 | 36 |
|
@@ -533,6 +535,160 @@ bool FITveto(T const& bc, DGCutparHolder const& diffCuts) |
533 | 535 | return false; |
534 | 536 | } |
535 | 537 |
|
| 538 | +inline void setBit(uint64_t w[4], int bit, bool val) |
| 539 | +{ |
| 540 | + if (!val) { |
| 541 | + return; |
| 542 | + } |
| 543 | + const int word = bit >> 6; |
| 544 | + const int offs = bit & 63; |
| 545 | + w[word] |= (static_cast<uint64_t>(1) << offs); |
| 546 | +} |
| 547 | + |
| 548 | +template <typename TFT0, typename TFV0A> |
| 549 | +inline void buildFT0FV0Words(TFT0 const& ft0, TFV0A const& fv0a, |
| 550 | + uint64_t thr1[4], uint64_t thr2[4], |
| 551 | + float thr1_FT0A = 25., float thr1_FT0C = 50., float thr1_FV0A = 50., |
| 552 | + float thr2_FT0A = 50., float thr2_FT0C = 100., float thr2_FV0A = 100.) |
| 553 | +{ |
| 554 | + thr1[0] = thr1[1] = thr1[2] = thr1[3] = 0ull; |
| 555 | + thr2[0] = thr2[1] = thr2[2] = thr2[3] = 0ull; |
| 556 | + |
| 557 | + constexpr int kFT0AOffset = 0; |
| 558 | + constexpr int kFT0COffset = 96; |
| 559 | + constexpr int kFV0Offset = 208; |
| 560 | + |
| 561 | + auto ampsA = ft0.amplitudeA(); |
| 562 | + const int nA = std::min<int>(ampsA.size(), 96); |
| 563 | + for (int i = 0; i < nA; ++i) { |
| 564 | + const auto a = ampsA[i]; |
| 565 | + setBit(thr1, kFT0AOffset + i, a >= thr1_FT0A); |
| 566 | + setBit(thr2, kFT0AOffset + i, a >= thr2_FT0A); |
| 567 | + } |
| 568 | + |
| 569 | + auto ampsC = ft0.amplitudeC(); |
| 570 | + const int nC = std::min<int>(ampsC.size(), 112); |
| 571 | + for (int i = 0; i < nC; ++i) { |
| 572 | + const auto a = ampsC[i]; |
| 573 | + setBit(thr1, kFT0COffset + i, a >= thr1_FT0C); |
| 574 | + setBit(thr2, kFT0COffset + i, a >= thr2_FT0C); |
| 575 | + } |
| 576 | + |
| 577 | + auto ampsV = fv0a.amplitude(); |
| 578 | + const int nV = std::min<int>(ampsV.size(), 48); |
| 579 | + for (int i = 0; i < nV; ++i) { |
| 580 | + const auto a = ampsV[i]; |
| 581 | + setBit(thr1, kFV0Offset + i, a >= thr1_FV0A); |
| 582 | + setBit(thr2, kFV0Offset + i, a >= thr2_FV0A); |
| 583 | + } |
| 584 | +} |
| 585 | + |
| 586 | +// ----------------------------------------------------------------------------- |
| 587 | +// return eta and phi of a given FIT channel based on the bitset |
| 588 | +// Bit layout contract: |
| 589 | +constexpr int kFT0Bits = 208; // FT0 total channels |
| 590 | +constexpr int kFV0Bits = 48; // FV0A channels |
| 591 | +constexpr int kTotalBits = 256; // 4*64 |
| 592 | + |
| 593 | +// FT0 side split |
| 594 | +constexpr int kFT0AChannels = 96; // FT0A channels are [0..95] |
| 595 | +constexpr int kFT0CChannels = 112; // FT0C channels are [96..207] |
| 596 | +static_assert(kFT0AChannels + kFT0CChannels == kFT0Bits); |
| 597 | + |
| 598 | +using Bits256 = std::array<uint64_t, 4>; |
| 599 | + |
| 600 | +inline Bits256 makeBits256(uint64_t w0, uint64_t w1, uint64_t w2, uint64_t w3) |
| 601 | +{ |
| 602 | + return {w0, w1, w2, w3}; |
| 603 | +} |
| 604 | + |
| 605 | +inline bool testBit(Bits256 const& w, int bit) |
| 606 | +{ |
| 607 | + if (bit < 0 || bit >= kTotalBits) { |
| 608 | + return false; |
| 609 | + } |
| 610 | + return (w[bit >> 6] >> (bit & 63)) & 1ULL; |
| 611 | +} |
| 612 | + |
| 613 | +struct FitBitRef { |
| 614 | + enum class Det : uint8_t { FT0, |
| 615 | + FV0, |
| 616 | + Unknown }; |
| 617 | + Det det = Det::Unknown; |
| 618 | + int ch = -1; // FT0: 0..207, FV0: 0..47 |
| 619 | + bool isC = false; // only meaningful for FT0 |
| 620 | +}; |
| 621 | + |
| 622 | +inline FitBitRef decodeFitBit(int bit) |
| 623 | +{ |
| 624 | + FitBitRef out; |
| 625 | + if (bit >= 0 && bit < kFT0Bits) { |
| 626 | + out.det = FitBitRef::Det::FT0; |
| 627 | + out.ch = bit; // FT0 channel id |
| 628 | + out.isC = (bit >= kFT0AChannels); // C side if in upper range |
| 629 | + return out; |
| 630 | + } |
| 631 | + if (bit >= kFT0Bits && bit < kTotalBits) { |
| 632 | + out.det = FitBitRef::Det::FV0; |
| 633 | + out.ch = bit - kFT0Bits; // FV0A channel id 0..47 |
| 634 | + return out; |
| 635 | + } |
| 636 | + return out; |
| 637 | +} |
| 638 | + |
| 639 | +template <typename FT0DetT, typename OffsetsT> |
| 640 | +inline double getPhiFT0_fromChannel(FT0DetT& ft0Det, int ft0Ch, OffsetsT const& offsetFT0, int i) |
| 641 | +{ |
| 642 | + ft0Det.calculateChannelCenter(); |
| 643 | + auto chPos = ft0Det.getChannelCenter(ft0Ch); |
| 644 | + |
| 645 | + const double x = chPos.X() + offsetFT0[i].getX(); |
| 646 | + const double y = chPos.Y() + offsetFT0[i].getY(); |
| 647 | + |
| 648 | + return RecoDecay::phi(x, y); |
| 649 | +} |
| 650 | + |
| 651 | +template <typename FT0DetT, typename OffsetsT> |
| 652 | +inline double getEtaFT0_fromChannel(FT0DetT& ft0Det, int ft0Ch, OffsetsT const& offsetFT0, int i) |
| 653 | +{ |
| 654 | + ft0Det.calculateChannelCenter(); |
| 655 | + auto chPos = ft0Det.getChannelCenter(ft0Ch); |
| 656 | + |
| 657 | + double x = chPos.X() + offsetFT0[i].getX(); |
| 658 | + double y = chPos.Y() + offsetFT0[i].getY(); |
| 659 | + double z = chPos.Z() + offsetFT0[i].getZ(); |
| 660 | + |
| 661 | + // If this channel belongs to FT0C, flip z (matches your original intent) |
| 662 | + const bool isC = (ft0Ch >= kFT0AChannels); |
| 663 | + if (isC) { |
| 664 | + z = -z; |
| 665 | + } |
| 666 | + |
| 667 | + const double r = std::sqrt(x * x + y * y); |
| 668 | + const double theta = std::atan2(r, z); |
| 669 | + return -std::log(std::tan(0.5 * theta)); |
| 670 | +} |
| 671 | + |
| 672 | +template <typename FT0DetT, typename OffsetsT> |
| 673 | +inline bool getPhiEtaFromFitBit(FT0DetT& ft0Det, |
| 674 | + int bit, |
| 675 | + OffsetsT const& offsetFT0, |
| 676 | + int iRunOffset, |
| 677 | + double& phi, |
| 678 | + double& eta) |
| 679 | +{ |
| 680 | + auto ref = decodeFitBit(bit); |
| 681 | + if (ref.det != FitBitRef::Det::FT0) { |
| 682 | + return false; |
| 683 | + } |
| 684 | + |
| 685 | + // FT0A: 0..95, FT0C: 96..207 |
| 686 | + const int ft0Ch = bit; |
| 687 | + phi = getPhiFT0_fromChannel(ft0Det, ft0Ch, offsetFT0, iRunOffset); |
| 688 | + eta = getEtaFT0_fromChannel(ft0Det, ft0Ch, offsetFT0, iRunOffset); |
| 689 | + return true; |
| 690 | +} |
| 691 | + |
536 | 692 | // ----------------------------------------------------------------------------- |
537 | 693 |
|
538 | 694 | template <typename T> |
|
0 commit comments