diff --git a/PWGJE/TableProducer/emcalCorrectionTask.cxx b/PWGJE/TableProducer/emcalCorrectionTask.cxx index 037f793a0d1..99806da08ac 100644 --- a/PWGJE/TableProducer/emcalCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalCorrectionTask.cxx @@ -140,6 +140,7 @@ struct EmcalCorrectionTask { Configurable applyGainCalibShift{"applyGainCalibShift", false, "Apply shift for cell gain calibration to use values before cell format change (Sept. 2023)"}; Configurable applySoftwareTriggerSelection{"applySoftwareTriggerSelection", false, "Apply software trigger selection"}; Configurable softwareTriggerSelection{"softwareTriggerSelection", "fGammaHighPtEMCAL,fGammaHighPtDCAL", "Default: fGammaHighPtEMCAL,fGammaHighPtDCAL"}; + Configurable storePerDFInfo{"storePerDFInfo", false, "store addition information per DF."}; // cross talk emulation configs EmcCrossTalkConf emcCrossTalkConf; @@ -198,6 +199,11 @@ struct EmcalCorrectionTask { static constexpr float TrackNotOnEMCal = -900.f; static constexpr int kMaxMatchesPerCluster = 20; // Maximum number of tracks to match per cluster + // cluster size + size_t nCluster = 0; + size_t nClusterAmb = 0; + size_t nCells = 0; + void init(InitContext const&) { LOG(debug) << "Start init!"; @@ -305,6 +311,7 @@ struct EmcalCorrectionTask { sigmaLongAxis{100, 0., 1.0, "#sigma^{2}_{long}"}, sigmaShortAxis{100, 0., 1.0, "#sigma^{2}_{short}"}, nCellAxis{60, -0.5, 59.5, "#it{n}_{cells}"}, + nClusterAxis{10001, -0.5, 10000.5, "#it{N}_{cluster}"}, energyDenseAxis = {7000, 0.f, 70.f, "#it{E}_{cell} (GeV)"}; o2::framework::AxisSpec axisDeltaEta{400, -0.2, 0.2, "#Delta#eta"}; o2::framework::AxisSpec axisDeltaPhi{400, -0.2, 0.2, "#Delta#varphi (rad)"}; @@ -379,6 +386,12 @@ struct EmcalCorrectionTask { mExtraTimeShiftRunRanges.emplace_back(536565, 536590); // Commisioning-LHC23r mExtraTimeShiftRunRanges.emplace_back(542280, 543854); // LHC23zv-LHC23zy mExtraTimeShiftRunRanges.emplace_back(559544, 559856); // PbPb 2024 + + if (storePerDFInfo.value) { + mHistManager.add("hNClusterDF", "hNClusterDF", O2HistType::kTH1D, {nClusterAxis}); + mHistManager.add("hNClusterAmbigousDF", "hNClusterAmbigousDF", O2HistType::kTH1D, {nClusterAxis}); + mHistManager.add("hNCellDF", "hNCellDF", O2HistType::kTH1D, {nClusterAxis}); + } } template @@ -404,6 +417,9 @@ struct EmcalCorrectionTask { int nCellsProcessed = 0; std::unordered_map numberCollsInBC; // Number of collisions mapped to the global BC index of all BCs std::unordered_map numberCellsInBC; // Number of cells mapped to the global BC index of all BCs to check whether EMCal was readout + nCluster = 0; + nClusterAmb = 0; + nCells = 0; for (const auto& bc : bcs) { LOG(debug) << "Next BC"; @@ -539,6 +555,11 @@ struct EmcalCorrectionTask { } // end of collision loop LOG(detail) << "Processed " << nBCsProcessed << " BCs with " << nCellsProcessed << " cells"; + if (storePerDFInfo) { + mHistManager.fill(HIST("hNClusterDF"), nCluster); + mHistManager.fill(HIST("hNClusterAmbigousDF"), nClusterAmb); + mHistManager.fill(HIST("hNCellDF"), nCells); + } } PROCESS_SWITCH(EmcalCorrectionTask, processFull, "run full analysis", true); @@ -551,6 +572,10 @@ struct EmcalCorrectionTask { int nCellsProcessed = 0; std::unordered_map numberCollsInBC; // Number of collisions mapped to the global BC index of all BCs std::unordered_map numberCellsInBC; // Number of cells mapped to the global BC index of all BCs to check whether EMCal was readout + + nCluster = 0; + nClusterAmb = 0; + nCells = 0; for (const auto& bc : bcs) { LOG(debug) << "Next BC"; @@ -690,6 +715,11 @@ struct EmcalCorrectionTask { } // end of collision loop LOG(detail) << "Processed " << nBCsProcessed << " BCs with " << nCellsProcessed << " cells"; + if (storePerDFInfo) { + mHistManager.fill(HIST("hNClusterDF"), nCluster); + mHistManager.fill(HIST("hNClusterAmbigousDF"), nClusterAmb); + mHistManager.fill(HIST("hNCellDF"), nCells); + } } PROCESS_SWITCH(EmcalCorrectionTask, processWithSecondaries, "run full analysis with secondary track matching", false); @@ -702,6 +732,11 @@ struct EmcalCorrectionTask { int nCellsProcessed = 0; std::unordered_map numberCollsInBC; // Number of collisions mapped to the global BC index of all BCs std::unordered_map numberCellsInBC; // Number of cells mapped to the global BC index of all BCs to check whether EMCal was readout + + nCluster = 0; + nClusterAmb = 0; + nCells = 0; + for (const auto& bc : bcs) { LOG(debug) << "Next BC"; // Convert aod::Calo to o2::emcal::Cell which can be used with the clusterizer. @@ -868,6 +903,11 @@ struct EmcalCorrectionTask { } // end of collision loop LOG(detail) << "Processed " << nBCsProcessed << " BCs with " << nCellsProcessed << " cells"; + if (storePerDFInfo) { + mHistManager.fill(HIST("hNClusterDF"), nCluster); + mHistManager.fill(HIST("hNClusterAmbigousDF"), nClusterAmb); + mHistManager.fill(HIST("hNCellDF"), nCells); + } } PROCESS_SWITCH(EmcalCorrectionTask, processMCFull, "run full analysis with MC info", false); @@ -880,6 +920,10 @@ struct EmcalCorrectionTask { int nCellsProcessed = 0; std::unordered_map numberCollsInBC; // Number of collisions mapped to the global BC index of all BCs std::unordered_map numberCellsInBC; // Number of cells mapped to the global BC index of all BCs to check whether EMCal was readout + + nCluster = 0; + nClusterAmb = 0; + nCells = 0; for (const auto& bc : bcs) { LOG(debug) << "Next BC"; // Convert aod::Calo to o2::emcal::Cell which can be used with the clusterizer. @@ -1049,6 +1093,11 @@ struct EmcalCorrectionTask { } // end of collision loop LOG(detail) << "Processed " << nBCsProcessed << " BCs with " << nCellsProcessed << " cells"; + if (storePerDFInfo) { + mHistManager.fill(HIST("hNClusterDF"), nCluster); + mHistManager.fill(HIST("hNClusterAmbigousDF"), nClusterAmb); + mHistManager.fill(HIST("hNCellDF"), nCells); + } } PROCESS_SWITCH(EmcalCorrectionTask, processMCWithSecondaries, "run full analysis with MC info", false); @@ -1059,6 +1108,10 @@ struct EmcalCorrectionTask { int nBCsProcessed = 0; int nCellsProcessed = 0; + nCluster = 0; + nClusterAmb = 0; + nCells = 0; + for (const auto& bc : bcs) { LOG(debug) << "Next BC"; // Convert aod::Calo to o2::emcal::Cell which can be used with the clusterizer. @@ -1165,6 +1218,11 @@ struct EmcalCorrectionTask { nBCsProcessed++; } // end of bc loop LOG(debug) << "Done with process BC."; + if (storePerDFInfo) { + mHistManager.fill(HIST("hNClusterDF"), nCluster); + mHistManager.fill(HIST("hNClusterAmbigousDF"), nClusterAmb); + mHistManager.fill(HIST("hNCellDF"), nCells); + } } PROCESS_SWITCH(EmcalCorrectionTask, processStandalone, "run stand alone analysis", false); @@ -1210,13 +1268,14 @@ struct EmcalCorrectionTask { void fillClusterTable(Collision const& col, math_utils::Point3D const& vertexPos, size_t iClusterizer, const gsl::span cellIndicesBC, MatchResult* indexMapPair = nullptr, const std::vector* trackGlobalIndex = nullptr, MatchResult* indexMapPairSecondaries = nullptr, const std::vector* secondariesGlobalIndex = nullptr) { // average number of cells per cluster, only used the reseve a reasonable amount for the clustercells table - const size_t nAvgNcells = 3; + // const size_t nAvgNcells = 3; // we found a collision, put the clusters into the none ambiguous table - clusters.reserve(mAnalysisClusters.size()); + clusters.reserve(nCluster + mAnalysisClusters.size()); if (!mClusterLabels.empty()) { - mcclusters.reserve(mClusterLabels.size()); + mcclusters.reserve(nCluster + mClusterLabels.size()); } - clustercells.reserve(mAnalysisClusters.size() * nAvgNcells); + // Since reserve triggers a fatal when its too small, it is not save for cells to use it unless we use a really large buffer... + // clustercells.reserve(mAnalysisClusters.size() * nAvgNcells); // get the clusterType once const auto clusterType = static_cast(mClusterDefinitions[iClusterizer]); @@ -1253,6 +1312,7 @@ struct EmcalCorrectionTask { cluster.getClusterTime(), cluster.getIsExotic(), cluster.getDistanceToBadChannel(), cluster.getNExMax(), clusterType); + ++nCluster; if (!mClusterLabels.empty()) { mcclusters(mClusterLabels[iCluster].getLabels(), mClusterLabels[iCluster].getEnergyFractions()); } @@ -1262,6 +1322,7 @@ struct EmcalCorrectionTask { LOG(debug) << "trying to find cell index " << cellindex << " in map"; if (cellIndicesBC[cellindex] >= 0) { clustercells(clusters.lastIndex(), cellIndicesBC[cellindex]); + ++nCells; } } // end of cells of cluser loop // fill histograms @@ -1305,13 +1366,13 @@ struct EmcalCorrectionTask { void fillAmbigousClusterTable(BC const& bc, size_t iClusterizer, const gsl::span cellIndicesBC, bool hasCollision) { // average number of cells per cluster, only used the reseve a reasonable amount for the clustercells table - const size_t nAvgNcells = 3; + // const size_t nAvgNcells = 3; int cellindex = -1; - clustersAmbiguous.reserve(mAnalysisClusters.size()); + clustersAmbiguous.reserve(mAnalysisClusters.size() + nClusterAmb); if (mClusterLabels.size() > 0) { - mcclustersAmbiguous.reserve(mClusterLabels.size()); + mcclustersAmbiguous.reserve(mClusterLabels.size() + nClusterAmb); } - clustercellsambiguous.reserve(mAnalysisClusters.size() * nAvgNcells); + // clustercellsambiguous.reserve(mAnalysisClusters.size() * nAvgNcells); unsigned int iCluster = 0; float energy = 0.f; for (const auto& cluster : mAnalysisClusters) { @@ -1343,6 +1404,7 @@ struct EmcalCorrectionTask { cluster.getM20(), cluster.getNCells(), cluster.getClusterTime(), cluster.getIsExotic(), cluster.getDistanceToBadChannel(), cluster.getNExMax(), static_cast(mClusterDefinitions.at(iClusterizer))); + ++nClusterAmb; if (mClusterLabels.size() > 0) { mcclustersAmbiguous(mClusterLabels[iCluster].getLabels(), mClusterLabels[iCluster].getEnergyFractions()); }