|
34 | 34 | #include <TPDGCode.h> |
35 | 35 |
|
36 | 36 | #include <algorithm> |
| 37 | +#include <cmath> |
37 | 38 | #include <string> |
| 39 | +#include <utility> |
38 | 40 | #include <vector> |
39 | 41 |
|
40 | 42 | using namespace o2; |
@@ -122,14 +124,49 @@ struct Cascqaanalysis { |
122 | 124 | SliceCache cache; |
123 | 125 |
|
124 | 126 | // Random number generator for event scaling |
125 | | - TRandom2* fRand = new TRandom2(); |
| 127 | + TRandom2 fRand; |
126 | 128 |
|
127 | 129 | // Struct to select on event type |
128 | 130 | typedef struct CollisionIndexAndType { |
129 | 131 | int64_t index; |
130 | 132 | uint8_t typeFlag; |
131 | 133 | } CollisionIndexAndType; |
132 | 134 |
|
| 135 | + template <typename TTrack> |
| 136 | + static int countITSHits(TTrack const& track) |
| 137 | + { |
| 138 | + int nHits = 0; |
| 139 | + for (unsigned int i = 0; i < 7; ++i) { |
| 140 | + if (track.itsClusterMap() & (1 << i)) { |
| 141 | + ++nHits; |
| 142 | + } |
| 143 | + } |
| 144 | + return nHits; |
| 145 | + } |
| 146 | + |
| 147 | + template <typename TCollision> |
| 148 | + static uint8_t buildRecoEventFlags(TCollision const& collision) |
| 149 | + { |
| 150 | + uint8_t evFlag = o2::aod::mycascades::EvFlags::EvINEL; |
| 151 | + if (collision.isInelGt0()) { |
| 152 | + evFlag |= o2::aod::mycascades::EvFlags::EvINELgt0; |
| 153 | + } |
| 154 | + if (collision.isInelGt1()) { |
| 155 | + evFlag |= o2::aod::mycascades::EvFlags::EvINELgt1; |
| 156 | + } |
| 157 | + return evFlag; |
| 158 | + } |
| 159 | + |
| 160 | + template <typename TCascade, typename TCollision> |
| 161 | + static std::pair<float, float> computeCascadeCtau(TCascade const& casc, TCollision const& collision) |
| 162 | + { |
| 163 | + const float decayLength = std::hypot(casc.x() - collision.posX(), casc.y() - collision.posY(), casc.z() - collision.posZ()); |
| 164 | + const float totalMomentum = std::hypot(casc.px(), casc.py(), casc.pz()); |
| 165 | + const float invMomentum = 1.f / (totalMomentum + 1.e-13f); |
| 166 | + return {o2::constants::physics::MassXiMinus * decayLength * invMomentum, |
| 167 | + o2::constants::physics::MassOmegaMinus * decayLength * invMomentum}; |
| 168 | + } |
| 169 | + |
133 | 170 | void init(InitContext const&) |
134 | 171 | { |
135 | 172 | TString hCandidateCounterLabels[4] = {"All candidates", "passed topo cuts", "has associated MC particle", "associated with Xi(Omega)"}; |
@@ -210,17 +247,13 @@ struct Cascqaanalysis { |
210 | 247 | auto bachelor = cascCand.template bachelor_as<TCascTracksTo>(); |
211 | 248 |
|
212 | 249 | // Basic set of selections |
213 | | - if (cascCand.cascradius() > cascradius && |
214 | | - cascCand.v0radius() > v0radius && |
215 | | - cascCand.casccosPA(pvx, pvy, pvz) > casccospa && |
216 | | - cascCand.v0cosPA(pvx, pvy, pvz) > v0cospa && |
217 | | - std::fabs(posdau.eta()) < etadau && |
218 | | - std::fabs(negdau.eta()) < etadau && |
219 | | - std::fabs(bachelor.eta()) < etadau) { |
220 | | - return true; |
221 | | - } else { |
222 | | - return false; |
223 | | - } |
| 250 | + return cascCand.cascradius() > cascradius && |
| 251 | + cascCand.v0radius() > v0radius && |
| 252 | + cascCand.casccosPA(pvx, pvy, pvz) > casccospa && |
| 253 | + cascCand.v0cosPA(pvx, pvy, pvz) > v0cospa && |
| 254 | + std::fabs(posdau.eta()) < etadau && |
| 255 | + std::fabs(negdau.eta()) < etadau && |
| 256 | + std::fabs(bachelor.eta()) < etadau; |
224 | 257 | } |
225 | 258 |
|
226 | 259 | template <typename TMcParticles> |
@@ -419,39 +452,16 @@ struct Cascqaanalysis { |
419 | 452 | registry.fill(HIST("hCandidateCounter"), 1.5); // passed topo cuts |
420 | 453 | nCandSel++; |
421 | 454 | // Fill table |
422 | | - if (fRand->Rndm() < lEventScale) { |
| 455 | + if (fRand.Rndm() < lEventScale) { |
423 | 456 | auto posdau = casc.posTrack_as<DauTracks>(); |
424 | 457 | auto negdau = casc.negTrack_as<DauTracks>(); |
425 | 458 | auto bachelor = casc.bachelor_as<DauTracks>(); |
426 | 459 |
|
427 | | - // ITS N hits |
428 | | - int posITSNhits = 0, negITSNhits = 0, bachITSNhits = 0; |
429 | | - for (unsigned int i = 0; i < 7; i++) { |
430 | | - if (posdau.itsClusterMap() & (1 << i)) { |
431 | | - posITSNhits++; |
432 | | - } |
433 | | - if (negdau.itsClusterMap() & (1 << i)) { |
434 | | - negITSNhits++; |
435 | | - } |
436 | | - if (bachelor.itsClusterMap() & (1 << i)) { |
437 | | - bachITSNhits++; |
438 | | - } |
439 | | - } |
440 | | - |
441 | | - uint8_t evFlag = 0; |
442 | | - evFlag |= o2::aod::mycascades::EvFlags::EvINEL; |
443 | | - if (collision.multNTracksPVeta1() > 0) { |
444 | | - evFlag |= o2::aod::mycascades::EvFlags::EvINELgt0; |
445 | | - } |
446 | | - if (collision.multNTracksPVeta1() > 1) { |
447 | | - evFlag |= o2::aod::mycascades::EvFlags::EvINELgt1; |
448 | | - } |
449 | | - |
450 | | - // c x tau |
451 | | - float cascpos = std::hypot(casc.x() - collision.posX(), casc.y() - collision.posY(), casc.z() - collision.posZ()); |
452 | | - float cascptotmom = std::hypot(casc.px(), casc.py(), casc.pz()); |
453 | | - float ctauXi = o2::constants::physics::MassXiMinus * cascpos / (cascptotmom + 1e-13); |
454 | | - float ctauOmega = o2::constants::physics::MassOmegaMinus * cascpos / (cascptotmom + 1e-13); |
| 460 | + const int posITSNhits = countITSHits(posdau); |
| 461 | + const int negITSNhits = countITSHits(negdau); |
| 462 | + const int bachITSNhits = countITSHits(bachelor); |
| 463 | + const uint8_t evFlag = buildRecoEventFlags(collision); |
| 464 | + const auto [ctauXi, ctauOmega] = computeCascadeCtau(casc, collision); |
455 | 465 |
|
456 | 466 | mycascades(collision.posZ(), |
457 | 467 | collision.centFT0M(), collision.centFV0A(), |
@@ -563,41 +573,17 @@ struct Cascqaanalysis { |
563 | 573 | genY = cascmc.y(); |
564 | 574 | } |
565 | 575 | } |
566 | | - if (fRand->Rndm() < lEventScale) { |
| 576 | + if (fRand.Rndm() < lEventScale) { |
567 | 577 | // Fill table |
568 | 578 | auto posdau = casc.posTrack_as<DauTracks>(); |
569 | 579 | auto negdau = casc.negTrack_as<DauTracks>(); |
570 | 580 | auto bachelor = casc.bachelor_as<DauTracks>(); |
571 | 581 |
|
572 | | - // ITS N hits |
573 | | - int posITSNhits = 0, negITSNhits = 0, bachITSNhits = 0; |
574 | | - for (unsigned int i = 0; i < 7; i++) { |
575 | | - if (posdau.itsClusterMap() & (1 << i)) { |
576 | | - posITSNhits++; |
577 | | - } |
578 | | - if (negdau.itsClusterMap() & (1 << i)) { |
579 | | - negITSNhits++; |
580 | | - } |
581 | | - if (bachelor.itsClusterMap() & (1 << i)) { |
582 | | - bachITSNhits++; |
583 | | - } |
584 | | - } |
585 | | - |
586 | | - // Event type flag |
587 | | - uint8_t evFlag = 0; |
588 | | - evFlag |= o2::aod::mycascades::EvFlags::EvINEL; |
589 | | - if (collision.multNTracksPVeta1() > 0) { |
590 | | - evFlag |= o2::aod::mycascades::EvFlags::EvINELgt0; |
591 | | - } |
592 | | - if (collision.multNTracksPVeta1() > 1) { |
593 | | - evFlag |= o2::aod::mycascades::EvFlags::EvINELgt1; |
594 | | - } |
595 | | - |
596 | | - // c x tau |
597 | | - float cascpos = std::hypot(casc.x() - collision.posX(), casc.y() - collision.posY(), casc.z() - collision.posZ()); |
598 | | - float cascptotmom = std::hypot(casc.px(), casc.py(), casc.pz()); |
599 | | - float ctauXi = o2::constants::physics::MassXiMinus * cascpos / (cascptotmom + 1e-13); |
600 | | - float ctauOmega = o2::constants::physics::MassOmegaMinus * cascpos / (cascptotmom + 1e-13); |
| 582 | + const int posITSNhits = countITSHits(posdau); |
| 583 | + const int negITSNhits = countITSHits(negdau); |
| 584 | + const int bachITSNhits = countITSHits(bachelor); |
| 585 | + const uint8_t evFlag = buildRecoEventFlags(collision); |
| 586 | + const auto [ctauXi, ctauOmega] = computeCascadeCtau(casc, collision); |
601 | 587 |
|
602 | 588 | mycascades(collision.posZ(), |
603 | 589 | mcCollision.centFT0M(), 0, // mcCollision.centFV0A() to be added |
|
0 commit comments