diff --git a/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx b/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx index f2aa922ec90..0f1d2103db5 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalPi0Qc.cxx @@ -28,6 +28,7 @@ #include "Common/CCDB/TriggerAliases.h" #include "Common/DataModel/EventSelection.h" +#include #include #include #include @@ -40,8 +41,11 @@ #include #include +#include +#include +#include // IWYU pragma: keep +#include #include -#include #include #include @@ -50,26 +54,23 @@ #include #include #include +#include #include #include #include #include -#include #include using namespace o2::framework; using namespace o2::framework::expressions; -using SelectedClusters = o2::soa::Filtered; using MyCollisions = o2::soa::Join; -using MyBCs = o2::soa::Join; -using SelectedAmbiguousClusters = o2::soa::Filtered; +using MyBCs = o2::soa::Join; struct Photon { - Photon(float eta_tmp, float phi_tmp, float energy_tmp, int clusid = 0) + Photon(float eta_tmp, float phi_tmp, float energy_tmp, int clusid = 0, uint8_t sm_tmp = 0) { eta = eta_tmp; phi = phi_tmp; - onDCal = (phi < 6 && phi > 4); energy = energy_tmp; theta = 2 * std::atan2(std::exp(-eta), 1); px = energy * std::sin(theta) * std::cos(phi); @@ -78,19 +79,22 @@ struct Photon { pt = std::sqrt(px * px + py * py); photon.SetPxPyPzE(px, py, pz, energy); id = clusid; + sm = sm_tmp; + onDCal = (phi < 6 && phi > 4); } - TLorentzVector photon; + ROOT::Math::PxPyPzEVector photon; float pt; float px; float py; float pz; float eta; float phi; - bool onDCal; // Checks whether photon is in phi region of the DCal, otherwise: EMCal float energy; float theta; int id; + uint8_t sm; + bool onDCal; // Checks whether photon is in phi region of the DCal, otherwise: EMCal }; struct Meson { @@ -101,11 +105,17 @@ struct Meson { } Photon pgamma1; Photon pgamma2; - TLorentzVector pMeson; + ROOT::Math::PxPyPzEVector pMeson; float getMass() const { return pMeson.M(); } float getPt() const { return pMeson.Pt(); } - float getOpeningAngle() const { return pgamma1.photon.Angle(pgamma2.photon.Vect()); } + float getOpeningAngle() const + { + float cosAngle = pgamma1.photon.Vect().Dot(pgamma2.photon.Vect()) / (pgamma1.photon.P() * pgamma2.photon.P()); + float angle = std::acos(std::clamp(cosAngle, -1.0f, 1.0f)); + return angle; + } + ROOT::Math::PxPyPzEVector getMathVector() const { return pMeson; } }; struct EventMixVec { @@ -156,8 +166,8 @@ struct EmcalPi0Qc { // define cluster filter. It selects only those clusters which are of the type // specified in the string mClusterDefinition,e.g. kV3Default, which is V3 clusterizer with default // clusterization parameters - o2::aod::EMCALClusterDefinition clusDef = o2::aod::emcalcluster::getClusterDefinitionFromString(mClusterDefinition.value); - Filter clusterDefinitionSelection = o2::aod::emcalcluster::definition == static_cast(clusDef); + // o2::aod::EMCALClusterDefinition clusDef = o2::aod::emcalcluster::getClusterDefinitionFromString(mClusterDefinition.value); + // Filter clusterDefinitionSelection = o2::aod::emcalcluster::definition == static_cast(clusDef); // define container for photons std::vector mPhotons; @@ -165,9 +175,17 @@ struct EmcalPi0Qc { // event mixing class EventMixVec evtMix; + o2::ccdb::CcdbApi ccdbApi; + int lastRunNumber = -1; // get the runnumber to obtain the SOR of the run to get t - SOR in (s) later + int64_t tsSOR = -1; + o2::aod::EMCALClusterDefinition clusDef = o2::aod::emcalcluster::kV3Default; + /// \brief Create output histograms and initialize geometry void init(InitContext const&) { + // init ccdb api + ccdbApi.init("https://alice-ccdb.cern.ch"); + // load geometry just in case we need it mGeometry = o2::emcal::Geometry::GetInstanceFromRunNumber(300000); @@ -175,20 +193,30 @@ struct EmcalPi0Qc { LOG(info) << "Creating histograms"; const AxisSpec bcAxis{3501, -0.5, 3500.5}; const AxisSpec energyAxis{makeClusterBinning(), "#it{E} (GeV)"}; + const AxisSpec collisionTimeAxis{1440, 0, 1440, "#it{t} - SOR (min)"}; + const AxisSpec clusterTimeAxis{1500, -600, 900, "t_{cl} (ns)"}; + const AxisSpec invmassAxis{invmassBinning, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; + const AxisSpec ptAxis{pTBinning, "#it{p}_{T} (GeV/#it{c})"}; + + if (doprocessCollision) { + mHistManager.add("events", "events;;#it{count}", HistType::kTH1F, {{7, 0.5, 7.5}}); + auto heventType = mHistManager.get(HIST("events")); + heventType->GetXaxis()->SetBinLabel(1, "All events"); + heventType->GetXaxis()->SetBinLabel(2, "sel8"); + heventType->GetXaxis()->SetBinLabel(3, "EMCal readout"); + heventType->GetXaxis()->SetBinLabel(4, "1+ Contributor"); + heventType->GetXaxis()->SetBinLabel(5, "z<10cm"); + heventType->GetXaxis()->SetBinLabel(6, "unique col"); + heventType->GetXaxis()->SetBinLabel(7, "EMCal cell>0"); + mHistManager.add("eventVertexZAll", "z-vertex of event (all events)", HistType::kTH1F, {{200, -20, 20}}); + mHistManager.add("eventVertexZSelected", "z-vertex of event (selected events)", HistType::kTH1F, {{200, -20, 20}}); + mHistManager.add("hEventPerTime", "number of events per time", HistType::kTH1F, {collisionTimeAxis}); + } - mHistManager.add("events", "events;;#it{count}", HistType::kTH1F, {{6, 0.5, 6.5}}); - auto heventType = mHistManager.get(HIST("events")); - heventType->GetXaxis()->SetBinLabel(1, "All events"); - heventType->GetXaxis()->SetBinLabel(2, "sel8 + readout"); - heventType->GetXaxis()->SetBinLabel(3, "1+ Contributor"); - heventType->GetXaxis()->SetBinLabel(4, "z<10cm"); - heventType->GetXaxis()->SetBinLabel(5, "unique col"); - heventType->GetXaxis()->SetBinLabel(6, "EMCAL cell>0"); - mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", HistType::kTH1F, {bcAxis}); - mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", HistType::kTH1F, {bcAxis}); - mHistManager.add("eventVertexZAll", "z-vertex of event (all events)", HistType::kTH1F, {{200, -20, 20}}); - mHistManager.add("eventVertexZSelected", "z-vertex of event (selected events)", HistType::kTH1F, {{200, -20, 20}}); - + if (doprocessAmbiguous) { + mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", HistType::kTH1F, {bcAxis}); + mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", HistType::kTH1F, {bcAxis}); + } // cluster properties for (const bool& iBeforeCuts : {false, true}) { const char* clusterDirectory = iBeforeCuts ? "ClustersBeforeCuts" : "ClustersAfterCuts"; @@ -204,17 +232,25 @@ struct EmcalPi0Qc { } // meson related histograms - mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPtMixedBackground", "invariant mass and pT of mixed background meson candidates", HistType::kTH2F, {invmassBinning, pTBinning}); + mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPtMixedBackground", "invariant mass and pT of mixed background meson candidates", HistType::kTH2F, {invmassAxis, ptAxis}); if (mSplitEMCalDCal) { - mHistManager.add("invMassVsPt_EMCal", "invariant mass and pT of meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPtBackground_EMCal", "invariant mass and pT of background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPtMixedBackground_EMCal", "invariant mass and pT of mixed background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPt_DCal", "invariant mass and pT of meson candidates with both clusters on DCal", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPtBackground_DCal", "invariant mass and pT of background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassBinning, pTBinning}); - mHistManager.add("invMassVsPtMixedBackground_DCal", "invariant mass and pT of mixed background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassBinning, pTBinning}); + mHistManager.add("invMassVsPt_EMCal", "invariant mass and pT of meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPtBackground_EMCal", "invariant mass and pT of background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPtMixedBackground_EMCal", "invariant mass and pT of mixed background meson candidates with both clusters on EMCal", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPt_DCal", "invariant mass and pT of meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPtBackground_DCal", "invariant mass and pT of background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}); + mHistManager.add("invMassVsPtMixedBackground_DCal", "invariant mass and pT of mixed background meson candidates with both clusters on DCal", HistType::kTH2F, {invmassAxis, ptAxis}); + } + + // add histograms per supermodule + for (int ism = 0; ism < 20; ++ism) { + mHistManager.add(Form("clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM%d", ism), Form("Cluster time vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {clusterTimeAxis, collisionTimeAxis}); + mHistManager.add(Form("clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM%d", ism), Form("Cluster number of cells vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {{50, 0, 50}, collisionTimeAxis}); + mHistManager.add(Form("clusterM02VsTimeStamp/clusterM02VsTimeStampSM%d", ism), Form("Cluster M02 vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {{400, 0, 5}, collisionTimeAxis}); + mHistManager.add(Form("mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM%d", ism), Form("invariant mass vs collision timestamp in Supermodule %d", ism), HistType::kTH2F, {invmassAxis, collisionTimeAxis}); } if (mVetoBCID->length()) { @@ -237,69 +273,258 @@ struct EmcalPi0Qc { mSelectBCIDs.push_back(bcid); } } + clusDef = o2::aod::emcalcluster::getClusterDefinitionFromString(mClusterDefinition.value); + LOG(info) << "mDoEventSel = " << mDoEventSel.value; + LOG(info) << "mRequireCaloReadout = " << mRequireCaloReadout.value; + LOG(info) << "mRequireEMCalCells = " << mRequireEMCalCells.value; + LOG(info) << "mSplitEMCalDCal = " << mSplitEMCalDCal.value; + } + + template + void supermoduleHistHelperPhoton(float time, float m02, int NCell, float timeSinceSOR) + { + static constexpr std::string_view ClusterTimeHistSM[20] = {"clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM0", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM1", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM2", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM3", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM4", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM5", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM6", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM7", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM8", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM9", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM10", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM11", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM12", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM13", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM14", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM15", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM16", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM17", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM18", "clusterTimeVsTimeStamp/clusterTimeVsTimeStampSM19"}; + static constexpr std::string_view ClusterNcellHistSM[20] = {"clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM0", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM1", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM2", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM3", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM4", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM5", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM6", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM7", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM8", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM9", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM10", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM11", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM12", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM13", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM14", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM15", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM16", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM17", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM18", "clusterNcellVsTimeStamp/clusterNCellVsTimeStampSM19"}; + static constexpr std::string_view ClusterM02HistSM[20] = {"clusterM02VsTimeStamp/clusterM02VsTimeStampSM0", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM1", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM2", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM3", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM4", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM5", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM6", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM7", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM8", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM9", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM10", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM11", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM12", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM13", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM14", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM15", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM16", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM17", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM18", "clusterM02VsTimeStamp/clusterM02VsTimeStampSM19"}; + mHistManager.fill(HIST(ClusterTimeHistSM[supermoduleID]), time, timeSinceSOR); + mHistManager.fill(HIST(ClusterNcellHistSM[supermoduleID]), NCell, timeSinceSOR); + mHistManager.fill(HIST(ClusterM02HistSM[supermoduleID]), m02, timeSinceSOR); + } + + template + void supermoduleHistHelperMeson(float minv, float timeSinceSOR) + { + static constexpr std::string_view MesonInvMassHistSM[20] = {"mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM0", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM1", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM2", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM3", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM4", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM5", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM6", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM7", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM8", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM9", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM10", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM11", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM12", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM13", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM14", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM15", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM16", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM17", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM18", "mesonInvMassVsTimeStamp/mesonInvMassVsTimeStampSM19"}; + mHistManager.fill(HIST(MesonInvMassHistSM[supermoduleID]), minv, timeSinceSOR); + } + + void fillSupermoduleHistogramsPhoton(int supermoduleID, float time, float m02, int NCell, float timeSinceSOR) + { + // workaround to have the histogram names per supermodule + // handled as constexpr + switch (supermoduleID) { + case 0: + supermoduleHistHelperPhoton<0>(time, m02, NCell, timeSinceSOR); + break; + case 1: + supermoduleHistHelperPhoton<1>(time, m02, NCell, timeSinceSOR); + break; + case 2: + supermoduleHistHelperPhoton<2>(time, m02, NCell, timeSinceSOR); + break; + case 3: + supermoduleHistHelperPhoton<3>(time, m02, NCell, timeSinceSOR); + break; + case 4: + supermoduleHistHelperPhoton<4>(time, m02, NCell, timeSinceSOR); + break; + case 5: + supermoduleHistHelperPhoton<5>(time, m02, NCell, timeSinceSOR); + break; + case 6: + supermoduleHistHelperPhoton<6>(time, m02, NCell, timeSinceSOR); + break; + case 7: + supermoduleHistHelperPhoton<7>(time, m02, NCell, timeSinceSOR); + break; + case 8: + supermoduleHistHelperPhoton<8>(time, m02, NCell, timeSinceSOR); + break; + case 9: + supermoduleHistHelperPhoton<9>(time, m02, NCell, timeSinceSOR); + break; + case 10: + supermoduleHistHelperPhoton<10>(time, m02, NCell, timeSinceSOR); + break; + case 11: + supermoduleHistHelperPhoton<11>(time, m02, NCell, timeSinceSOR); + break; + case 12: + supermoduleHistHelperPhoton<12>(time, m02, NCell, timeSinceSOR); + break; + case 13: + supermoduleHistHelperPhoton<13>(time, m02, NCell, timeSinceSOR); + break; + case 14: + supermoduleHistHelperPhoton<14>(time, m02, NCell, timeSinceSOR); + break; + case 15: + supermoduleHistHelperPhoton<15>(time, m02, NCell, timeSinceSOR); + break; + case 16: + supermoduleHistHelperPhoton<16>(time, m02, NCell, timeSinceSOR); + break; + case 17: + supermoduleHistHelperPhoton<17>(time, m02, NCell, timeSinceSOR); + break; + case 18: + supermoduleHistHelperPhoton<18>(time, m02, NCell, timeSinceSOR); + break; + case 19: + supermoduleHistHelperPhoton<19>(time, m02, NCell, timeSinceSOR); + break; + default: + break; + } + } + + void fillSupermoduleHistogramsMeson(int supermoduleID, float minv, float timeSinceSOR) + { + // workaround to have the histogram names per supermodule + // handled as constexpr + switch (supermoduleID) { + case 0: + supermoduleHistHelperMeson<0>(minv, timeSinceSOR); + break; + case 1: + supermoduleHistHelperMeson<1>(minv, timeSinceSOR); + break; + case 2: + supermoduleHistHelperMeson<2>(minv, timeSinceSOR); + break; + case 3: + supermoduleHistHelperMeson<3>(minv, timeSinceSOR); + break; + case 4: + supermoduleHistHelperMeson<4>(minv, timeSinceSOR); + break; + case 5: + supermoduleHistHelperMeson<5>(minv, timeSinceSOR); + break; + case 6: + supermoduleHistHelperMeson<6>(minv, timeSinceSOR); + break; + case 7: + supermoduleHistHelperMeson<7>(minv, timeSinceSOR); + break; + case 8: + supermoduleHistHelperMeson<8>(minv, timeSinceSOR); + break; + case 9: + supermoduleHistHelperMeson<9>(minv, timeSinceSOR); + break; + case 10: + supermoduleHistHelperMeson<10>(minv, timeSinceSOR); + break; + case 11: + supermoduleHistHelperMeson<11>(minv, timeSinceSOR); + break; + case 12: + supermoduleHistHelperMeson<12>(minv, timeSinceSOR); + break; + case 13: + supermoduleHistHelperMeson<13>(minv, timeSinceSOR); + break; + case 14: + supermoduleHistHelperMeson<14>(minv, timeSinceSOR); + break; + case 15: + supermoduleHistHelperMeson<15>(minv, timeSinceSOR); + break; + case 16: + supermoduleHistHelperMeson<16>(minv, timeSinceSOR); + break; + case 17: + supermoduleHistHelperMeson<17>(minv, timeSinceSOR); + break; + case 18: + supermoduleHistHelperMeson<18>(minv, timeSinceSOR); + break; + case 19: + supermoduleHistHelperMeson<19>(minv, timeSinceSOR); + break; + default: + break; + } } - PresliceUnsorted perCollision = o2::aod::emcalcluster::collisionId; + PresliceUnsortedOptional perCollision = o2::aod::emcalcluster::collisionId; + PresliceOptional perCluster = o2::aod::emcalclustercell::emcalclusterId; /// \brief Process EMCAL clusters that are matched to a collisions - void processCollision(MyBCs const&, MyCollisions const& collisions, SelectedClusters const& clusters, o2::soa::Filtered const& cells) + void processCollision(MyBCs const& bcs, MyCollisions const& collisions, o2::aod::EMCALClusters const& clusters, o2::soa::Filtered const& cells, o2::aod::EMCALClusterCells const& clusterCells) { + + auto cellIter = cells.begin(); + auto bcIter = bcs.begin(); + int runNumber = bcIter.runNumber(); std::unordered_map cellGlobalBCs; // Build map of number of cells for corrected BCs using global BCs // used later in the determination whether a BC has EMC cell content (for speed reason) for (const auto& cell : cells) { - auto globalbcid = cell.bc_as().globalBC(); - auto found = cellGlobalBCs.find(globalbcid); - if (found != cellGlobalBCs.end()) { - found->second++; - } else { - cellGlobalBCs.insert(std::pair(globalbcid, 1)); - } + cellGlobalBCs[cell.bc_as().globalBC()]++; } for (const auto& collision : collisions) { mHistManager.fill(HIST("events"), 1); // Fill "All events" bin of event histogram - if (mDoEventSel && (!collision.sel8() || (mRequireCaloReadout && !collision.alias_bit(kTVXinEMC)))) { // Check sel8 and whether EMC was read out + if (mDoEventSel.value && (!collision.sel8())) { // Check sel8 + continue; + } + mHistManager.fill(HIST("events"), 2); // Fill sel8 + if (mRequireCaloReadout.value && !collision.alias_bit(kTVXinEMC)) { // Check whether EMC was read out continue; } - mHistManager.fill(HIST("events"), 2); // Fill sel8 + readout + mHistManager.fill(HIST("events"), 3); // Fill readout - if (mDoEventSel && collision.numContrib() < 0.5) { // Skip collisions without contributors + if (mDoEventSel.value && collision.numContrib() < 0.5) { // Skip collisions without contributors continue; } - mHistManager.fill(HIST("events"), 3); // Fill >1 vtx contr. bin of event histogram + mHistManager.fill(HIST("events"), 4); // Fill >1 vtx contr. bin of event histogram mHistManager.fill(HIST("eventVertexZAll"), collision.posZ()); if (mVertexCut > 0 && std::abs(collision.posZ()) > mVertexCut) { continue; } - mHistManager.fill(HIST("events"), 4); // Fill z-Vertex selected bin of event histogram + mHistManager.fill(HIST("events"), 5); // Fill z-Vertex selected bin of event histogram mHistManager.fill(HIST("eventVertexZSelected"), collision.posZ()); - if (mDoEventSel && collision.ambiguous()) { // Skip ambiguous collisions (those that are in BCs including multiple collisions) + if (mDoEventSel.value && collision.ambiguous()) { // Skip ambiguous collisions (those that are in BCs including multiple collisions) continue; } - mHistManager.fill(HIST("events"), 5); // Fill "One collision in BC" bin of event histogram + mHistManager.fill(HIST("events"), 6); // Fill "One collision in BC" bin of event histogram - if (mDoEventSel) { + if (mDoEventSel.value) { auto found = cellGlobalBCs.find(collision.foundBC_as().globalBC()); - if (mRequireEMCalCells && (found == cellGlobalBCs.end() || found->second == 0)) { // Skip collisions without any readout EMCal cells + if (mRequireEMCalCells.value && (found == cellGlobalBCs.end() || found->second == 0)) { // Skip collisions without any readout EMCal cells continue; } } - mHistManager.fill(HIST("events"), 6); // Fill at least one non0 cell in EMCal of event histogram (Selected) + mHistManager.fill(HIST("events"), 7); // Fill at least one non0 cell in EMCal of event histogram (Selected) - auto clustersPerColl = clusters.sliceBy(perCollision, collision.collisionId()); - processClusters(clustersPerColl); - processMesons(); + // Get BC and run number + int64_t foundBCId = collision.foundBCId(); + if (foundBCId >= 0) { + bcIter.setCursor(foundBCId); + } + runNumber = bcIter.runNumber(); + + // Fetch SOR only when run changes + if (runNumber != lastRunNumber) { + std::map headers, metadata; + headers = ccdbApi.retrieveHeaders(Form("RCT/Info/RunInformation/%i", runNumber), metadata, -1); + tsSOR = atol(headers["SOR"].c_str()); + // LOGP(info, "Run {} | SOR = {} ms", runNumber, tsSOR); + lastRunNumber = runNumber; + } + + // Time since SOR in minutes (bc.timestamp() is in ms) + float timeSinceSORMin = (bcIter.timestamp() - tsSOR) / 1000.0f / 60.f; + mHistManager.fill(HIST("hEventPerTime"), timeSinceSORMin); + + auto clustersPerColl = clusters.sliceBy(perCollision, collision.globalIndex()); + if (clustersPerColl.size() == 0) { + continue; + } + processClusters(clustersPerColl, clusterCells, cellIter, timeSinceSORMin); + processMesons(timeSinceSORMin); } } PROCESS_SWITCH(EmcalPi0Qc, processCollision, "Process clusters from collision", false); /// \brief Process EMCAL clusters that are not matched to a collision /// This is not needed for most users - void processAmbiguous(o2::aod::BCs::iterator const& bc, SelectedAmbiguousClusters const& clusters) + void processAmbiguous(o2::aod::BCs::iterator const& bc, o2::aod::EMCALAmbiguousClusters const& clusters) { LOG(debug) << "processAmbiguous"; // TODO: remove this loop and put it in separate process function that only takes care of ambiguous clusters @@ -321,42 +546,39 @@ struct EmcalPi0Qc { PROCESS_SWITCH(EmcalPi0Qc, processAmbiguous, "Process Ambiguous clusters", false); /// \brief Process EMCAL clusters that are matched to a collisions - template - void processClusters(Clusters const& clusters) + template + void processClusters(Clusters const& clusters, o2::aod::EMCALClusterCells const& clusterCells, Cell& cellIter, float timeSinceSOR = 0.f) { LOG(debug) << "processClusters"; // clear photon vector mPhotons.clear(); - int globalCollID = -1000; - - // loop over all clusters from accepted collision - // auto eventClusters = clusters.select(o2::aod::emcalcluster::bcId == theCollision.bc().globalBC()); for (const auto& cluster : clusters) { - - // o2::InteractionRecord eventIR; - auto collID = cluster.collisionId(); - if (globalCollID == -1000) - globalCollID = collID; - - if (globalCollID != collID) { - LOG(info) << "Something went wrong with the collision ID"; + if (static_cast(clusDef) != cluster.definition()) { + continue; } + auto cellsPerCluster = clusterCells.sliceBy(perCluster, cluster.globalIndex()); + auto cellsPerClusterIter = cellsPerCluster.begin(); + cellIter.setCursor(cellsPerClusterIter.caloId()); + auto [supermodule, module, phiInModule, etaInModule] = mGeometry->GetCellIndex(cellIter.cellNumber()); + fillClusterQAHistos(cluster); - if (clusterRejectedByCut(cluster)) + if (clusterRejectedByCut(cluster)) { continue; + } fillClusterQAHistos(cluster); // put clusters in photon vector - mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy(), cluster.id())); + fillSupermoduleHistogramsPhoton(supermodule, cluster.time(), cluster.m02(), cluster.nCells(), timeSinceSOR); + mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy(), cluster.id(), supermodule)); } } /// \brief Process EMCAL clusters that are not matched to a collisions - template + template void processAmbiguousClusters(Clusters const& clusters) { LOG(debug) << "processClusters"; @@ -379,33 +601,33 @@ struct EmcalPi0Qc { } /// \brief Fills the standard QA histograms for a given cluster - template + template void fillClusterQAHistos(Cluster const& cluster) { // In this implementation the cluster properties are directly loaded from the flat table, // in the future one should consider using the AnalysisCluster object to work with after loading. - static constexpr std::string_view kClusterQAHistEnergy[2] = {"ClustersBeforeCuts/clusterE", "ClustersAfterCuts/clusterE"}; - static constexpr std::string_view kClusterQAHistEnergySimpleBinning[2] = {"ClustersBeforeCuts/clusterE_SimpleBinning", "ClustersAfterCuts/clusterE_SimpleBinning"}; - static constexpr std::string_view kClusterQAHistTime[2] = {"ClustersBeforeCuts/clusterTime", "ClustersAfterCuts/clusterTime"}; - static constexpr std::string_view kClusterQAHistEtaPhi[2] = {"ClustersBeforeCuts/clusterEtaPhi", "ClustersAfterCuts/clusterEtaPhi"}; - static constexpr std::string_view kClusterQAHistM02[2] = {"ClustersBeforeCuts/clusterM02", "ClustersAfterCuts/clusterM02"}; - static constexpr std::string_view kClusterQAHistM20[2] = {"ClustersBeforeCuts/clusterM20", "ClustersAfterCuts/clusterM20"}; - static constexpr std::string_view kClusterQAHistNLM[2] = {"ClustersBeforeCuts/clusterNLM", "ClustersAfterCuts/clusterNLM"}; - static constexpr std::string_view kClusterQAHistNCells[2] = {"ClustersBeforeCuts/clusterNCells", "ClustersAfterCuts/clusterNCells"}; - static constexpr std::string_view kClusterQAHistDistanceToBadChannel[2] = {"ClustersBeforeCuts/clusterDistanceToBadChannel", "ClustersAfterCuts/clusterDistanceToBadChannel"}; - mHistManager.fill(HIST(kClusterQAHistEnergy[BeforeCuts]), cluster.energy()); - mHistManager.fill(HIST(kClusterQAHistEnergySimpleBinning[BeforeCuts]), cluster.energy()); - mHistManager.fill(HIST(kClusterQAHistTime[BeforeCuts]), cluster.time()); - mHistManager.fill(HIST(kClusterQAHistEtaPhi[BeforeCuts]), cluster.eta(), cluster.phi()); - mHistManager.fill(HIST(kClusterQAHistM02[BeforeCuts]), cluster.m02()); - mHistManager.fill(HIST(kClusterQAHistM20[BeforeCuts]), cluster.m20()); - mHistManager.fill(HIST(kClusterQAHistNLM[BeforeCuts]), cluster.nlm()); - mHistManager.fill(HIST(kClusterQAHistNCells[BeforeCuts]), cluster.nCells()); - mHistManager.fill(HIST(kClusterQAHistDistanceToBadChannel[BeforeCuts]), cluster.distanceToBadChannel()); + static constexpr std::string_view ClusterQAHistEnergy[2] = {"ClustersBeforeCuts/clusterE", "ClustersAfterCuts/clusterE"}; + static constexpr std::string_view ClusterQAHistEnergySimpleBinning[2] = {"ClustersBeforeCuts/clusterE_SimpleBinning", "ClustersAfterCuts/clusterE_SimpleBinning"}; + static constexpr std::string_view ClusterQAHistTime[2] = {"ClustersBeforeCuts/clusterTime", "ClustersAfterCuts/clusterTime"}; + static constexpr std::string_view ClusterQAHistEtaPhi[2] = {"ClustersBeforeCuts/clusterEtaPhi", "ClustersAfterCuts/clusterEtaPhi"}; + static constexpr std::string_view ClusterQAHistM02[2] = {"ClustersBeforeCuts/clusterM02", "ClustersAfterCuts/clusterM02"}; + static constexpr std::string_view ClusterQAHistM20[2] = {"ClustersBeforeCuts/clusterM20", "ClustersAfterCuts/clusterM20"}; + static constexpr std::string_view ClusterQAHistNLM[2] = {"ClustersBeforeCuts/clusterNLM", "ClustersAfterCuts/clusterNLM"}; + static constexpr std::string_view ClusterQAHistNCells[2] = {"ClustersBeforeCuts/clusterNCells", "ClustersAfterCuts/clusterNCells"}; + static constexpr std::string_view ClusterQAHistDistanceToBadChannel[2] = {"ClustersBeforeCuts/clusterDistanceToBadChannel", "ClustersAfterCuts/clusterDistanceToBadChannel"}; + mHistManager.fill(HIST(ClusterQAHistEnergy[BeforeCuts]), cluster.energy()); + mHistManager.fill(HIST(ClusterQAHistEnergySimpleBinning[BeforeCuts]), cluster.energy()); + mHistManager.fill(HIST(ClusterQAHistTime[BeforeCuts]), cluster.time()); + mHistManager.fill(HIST(ClusterQAHistEtaPhi[BeforeCuts]), cluster.eta(), cluster.phi()); + mHistManager.fill(HIST(ClusterQAHistM02[BeforeCuts]), cluster.m02()); + mHistManager.fill(HIST(ClusterQAHistM20[BeforeCuts]), cluster.m20()); + mHistManager.fill(HIST(ClusterQAHistNLM[BeforeCuts]), cluster.nlm()); + mHistManager.fill(HIST(ClusterQAHistNCells[BeforeCuts]), cluster.nCells()); + mHistManager.fill(HIST(ClusterQAHistDistanceToBadChannel[BeforeCuts]), cluster.distanceToBadChannel()); } /// \brief Return a boolean that states, whether a cluster should be rejected by the applied cluster cuts - template + template bool clusterRejectedByCut(Cluster const& cluster) { // apply basic cluster cuts @@ -432,7 +654,7 @@ struct EmcalPi0Qc { } /// \brief Process meson candidates, calculate invariant mass and pT and fill histograms - void processMesons() + void processMesons(float timeSinceSOR = 0.f) { LOG(debug) << "processMesons " << mPhotons.size(); @@ -449,6 +671,12 @@ struct EmcalPi0Qc { if (meson.getOpeningAngle() > mMinOpenAngleCut) { mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt()); + uint8_t sm1 = mPhotons[ig1].sm; + uint8_t sm2 = mPhotons[ig2].sm; + if (sm1 == sm2) { + fillSupermoduleHistogramsMeson(sm1, meson.getMass(), timeSinceSOR); + } + if (mSplitEMCalDCal) { if (!mPhotons[ig1].onDCal && !mPhotons[ig2].onDCal) { mHistManager.fill(HIST("invMassVsPt_EMCal"), meson.getMass(), meson.getPt()); @@ -476,24 +704,28 @@ struct EmcalPi0Qc { } const double rotationAngle = o2::constants::math::PIHalf; // 0.78539816339; // rotaion angle 90° - TLorentzVector lvRotationPhoton1; // photon candidates which get rotated - TLorentzVector lvRotationPhoton2; // photon candidates which get rotated - TVector3 lvRotationPion; // rotation axis + ROOT::Math::PxPyPzEVector lvRotationPhoton1; // photon candidates which get rotated + ROOT::Math::PxPyPzEVector lvRotationPhoton2; // photon candidates which get rotated + ROOT::Math::PxPyPzEVector lvRotationPion; // rotation axis for (unsigned int ig3 = 0; ig3 < mPhotons.size(); ++ig3) { // continue if photons are identical if (ig3 == ig1 || ig3 == ig2) { continue; } - // calculate rotation axis - lvRotationPion = (meson.pMeson).Vect(); // initialize photons for rotation lvRotationPhoton1.SetPxPyPzE(mPhotons[ig1].px, mPhotons[ig1].py, mPhotons[ig1].pz, mPhotons[ig1].energy); lvRotationPhoton2.SetPxPyPzE(mPhotons[ig2].px, mPhotons[ig2].py, mPhotons[ig2].pz, mPhotons[ig2].energy); + lvRotationPion = meson.getMathVector(); + + // calculate rotation axis and matrix + lvRotationPion = lvRotationPhoton1 + lvRotationPhoton2; + ROOT::Math::AxisAngle rotationAxis(lvRotationPion.Vect(), rotationAngle); + ROOT::Math::Rotation3D rotationMatrix(rotationAxis); // rotate photons around rotation axis - lvRotationPhoton1.Rotate(rotationAngle, lvRotationPion); - lvRotationPhoton2.Rotate(rotationAngle, lvRotationPion); + lvRotationPhoton1 = rotationMatrix * lvRotationPhoton1; + lvRotationPhoton2 = rotationMatrix * lvRotationPhoton2; // initialize Photon objects for rotated photons Photon rotPhoton1(lvRotationPhoton1.Eta(), lvRotationPhoton1.Phi(), lvRotationPhoton1.E(), mPhotons[ig1].id);