From 2529cd6055e5af1b9a946628ef0a0334c5930d28 Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Thu, 7 Jul 2022 12:28:15 +0200 Subject: [PATCH 1/4] [EMCAL-734] adapted cluster definition --- PWGJE/DataModel/EMCALClusterDefinition.h | 8 ++++--- PWGJE/DataModel/EMCALClusters.h | 12 +++++------ PWGJE/TableProducer/emcalCorrectionTask.cxx | 5 ++++- PWGJE/Tasks/emcclustermonitor.cxx | 24 +++++++++------------ 4 files changed, 25 insertions(+), 24 deletions(-) diff --git a/PWGJE/DataModel/EMCALClusterDefinition.h b/PWGJE/DataModel/EMCALClusterDefinition.h index fee26cd8071..df2a1cdd012 100644 --- a/PWGJE/DataModel/EMCALClusterDefinition.h +++ b/PWGJE/DataModel/EMCALClusterDefinition.h @@ -39,12 +39,13 @@ struct EMCALClusterDefinition { double minCellEnergy = 0.05; // minimum cell energy (GeV) double timeMin = -10000; // minimum time (ns) double timeMax = 10000; // maximum time (ns) - double gradientCut = 0.03; // gradient cut + bool doGradientCut = true; // apply gradient cut if true + double gradientCut = -1; // gradient cut // default constructor EMCALClusterDefinition() = default; // constructor - EMCALClusterDefinition(ClusterAlgorithm_t pAlgorithm, int pStorageID, int pSelectedCellType, std::string pName, double pSeedEnergy, double pMinCellEnergy, double pTimeMin, double pTimeMax, double pGradientCut) + EMCALClusterDefinition(ClusterAlgorithm_t pAlgorithm, int pStorageID, int pSelectedCellType, std::string pName, double pSeedEnergy, double pMinCellEnergy, double pTimeMin, double pTimeMax, bool pDoGradientCut, double pGradientCut) { algorithm = pAlgorithm; storageID = pStorageID; @@ -54,13 +55,14 @@ struct EMCALClusterDefinition { minCellEnergy = pMinCellEnergy; timeMin = pTimeMin; timeMax = pTimeMax; + doGradientCut = pDoGradientCut; gradientCut = pGradientCut; } // implement comparison operators for int std::string and ClusterAlgorithm_t bool operator==(const EMCALClusterDefinition& rhs) const { - return (algorithm == rhs.algorithm && storageID == rhs.storageID && name == rhs.name && seedEnergy == rhs.seedEnergy && minCellEnergy == rhs.minCellEnergy && timeMin == rhs.timeMin && timeMax == rhs.timeMax && gradientCut == rhs.gradientCut); + return (algorithm == rhs.algorithm && storageID == rhs.storageID && name == rhs.name && seedEnergy == rhs.seedEnergy && minCellEnergy == rhs.minCellEnergy && timeMin == rhs.timeMin && timeMax == rhs.timeMax && gradientCut == rhs.gradientCut && doGradientCut == rhs.doGradientCut); } bool operator!=(const EMCALClusterDefinition& rhs) const { diff --git a/PWGJE/DataModel/EMCALClusters.h b/PWGJE/DataModel/EMCALClusters.h index 1969258ece2..8dc8f5d18c7 100644 --- a/PWGJE/DataModel/EMCALClusters.h +++ b/PWGJE/DataModel/EMCALClusters.h @@ -27,12 +27,12 @@ namespace emcalcluster // define global cluster definitions // the V1 algorithm is not yet implemented, but the V3 algorithm is // New definitions should be added here! -const EMCALClusterDefinition kV1Default(ClusterAlgorithm_t::kV1, 0, 1, "kV1Default", 0.5, 0.1, -10000, 10000, 0.03); -const EMCALClusterDefinition kV1Variation1(ClusterAlgorithm_t::kV3, 1, 1, "kV1Variation1", 0.3, 0.1, -10000, 10000, 0.03); -const EMCALClusterDefinition kV1Variation2(ClusterAlgorithm_t::kV3, 2, 1, "kV1Variation2", 0.2, 0.1, -10000, 10000, 0.03); -const EMCALClusterDefinition kV3Default(ClusterAlgorithm_t::kV3, 10, 1, "kV3Default", 0.5, 0.1, -10000, 10000, 0.03); -const EMCALClusterDefinition kV3Variation1(ClusterAlgorithm_t::kV3, 11, 1, "kV3Variation1", 0.3, 0.1, -10000, 10000, 0.03); -const EMCALClusterDefinition kV3Variation2(ClusterAlgorithm_t::kV3, 12, 1, "kV3Variation2", 0.2, 0.1, -10000, 10000, 0.03); +const EMCALClusterDefinition kV1Default(ClusterAlgorithm_t::kV1, 0, 1, "kV1Default", 0.5, 0.1, -10000, 10000, true, 0.03); //dummy +const EMCALClusterDefinition kV1Variation1(ClusterAlgorithm_t::kV3, 1, 1, "kV1Variation1", 0.3, 0.1, -10000, 10000, true, 0.03); //dummy +const EMCALClusterDefinition kV1Variation2(ClusterAlgorithm_t::kV3, 2, 1, "kV1Variation2", 0.2, 0.1, -10000, 10000, true, 0.03); //dummy +const EMCALClusterDefinition kV3Default(ClusterAlgorithm_t::kV3, 10, 1, "kV3Default", 0.5, 0.1, -10000, 10000, true, 0.03); +const EMCALClusterDefinition kV3Variation1(ClusterAlgorithm_t::kV3, 11, 1, "kV3Variation1", 0.5, 0.1, -10000, 10000, true, 0.); +const EMCALClusterDefinition kV3Variation2(ClusterAlgorithm_t::kV3, 12, 1, "kV3Variation2", 0.5, 0.1, -10000, 10000, false, 0.); /// \brief function returns EMCALClusterDefinition for the given name /// \param name name of the cluster definition diff --git a/PWGJE/TableProducer/emcalCorrectionTask.cxx b/PWGJE/TableProducer/emcalCorrectionTask.cxx index adb9515eddc..f82060f53d9 100644 --- a/PWGJE/TableProducer/emcalCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalCorrectionTask.cxx @@ -100,7 +100,7 @@ struct EmcalCorrectionTask { } } for (auto& clusterDefinition : mClusterDefinitions) { - mClusterizers.emplace_back(std::make_unique>(1, clusterDefinition.timeMin, clusterDefinition.timeMax, clusterDefinition.gradientCut, true, clusterDefinition.seedEnergy, clusterDefinition.minCellEnergy)); + mClusterizers.emplace_back(std::make_unique>(1E9, clusterDefinition.timeMin, clusterDefinition.timeMax, clusterDefinition.gradientCut, clusterDefinition.doGradientCut, clusterDefinition.seedEnergy, clusterDefinition.minCellEnergy)); mClusterFactories.emplace_back(std::make_unique>()); LOG(info) << "Cluster definition initialized: " << clusterDefinition.toString(); LOG(info) << "timeMin: " << clusterDefinition.timeMin; @@ -108,6 +108,7 @@ struct EmcalCorrectionTask { LOG(info) << "gradientCut: " << clusterDefinition.gradientCut; LOG(info) << "seedEnergy: " << clusterDefinition.seedEnergy; LOG(info) << "minCellEnergy: " << clusterDefinition.minCellEnergy; + LOG(info) << "storageID" << clusterDefinition.storageID; } for (auto& clusterizer : mClusterizers) { clusterizer->setGeometry(geometry); @@ -153,6 +154,7 @@ struct EmcalCorrectionTask { cell.time(), o2::emcal::intToChannelType(cell.cellType()))); } + LOG(debug) << "Number of cells (CF): " << mEmcalCells.size(); // Cell QA // For convenience, use the clusterizer stored geometry to get the eta-phi @@ -197,6 +199,7 @@ struct EmcalCorrectionTask { for (int icl = 0; icl < mClusterFactories.at(i)->getNumberOfClusters(); icl++) { auto analysisCluster = mClusterFactories.at(i)->buildCluster(icl); mAnalysisClusters.emplace_back(analysisCluster); + LOG(debug) << "Cluster " << icl << ": E: " << analysisCluster.E() << ", NCells " << analysisCluster.getNCells(); } LOG(debug) << "Converted to analysis clusters."; diff --git a/PWGJE/Tasks/emcclustermonitor.cxx b/PWGJE/Tasks/emcclustermonitor.cxx index 3c02bf5802c..03d3f7e5f6f 100644 --- a/PWGJE/Tasks/emcclustermonitor.cxx +++ b/PWGJE/Tasks/emcclustermonitor.cxx @@ -65,16 +65,11 @@ struct ClusterMonitor { Configurable mVetoBCID{"vetoBCID", "", "BC ID(s) to be excluded, this should be used as an alternative to the event selection"}; Configurable mSelectBCID{"selectBCID", "all", "BC ID(s) to be included, this should be used as an alternative to the event selection"}; Configurable mVertexCut{"vertexCut", -1, "apply z-vertex cut with value in cm"}; - Configurable mClusterDefinition{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"}; + Configurable mClusterDefinition{"clusterDefinition", 10, "cluster definition to be selected, e.g. 10=kV3Default"}; + std::vector mVetoBCIDs; std::vector mSelectBCIDs; - // 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); - /// \brief Create output histograms and initialize geometry void init(InitContext const&) { @@ -105,7 +100,7 @@ struct ClusterMonitor { // cluster properties (matched clusters) mHistManager.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis}); - mHistManager.add("clusterE_SimpleBinning", "Energy of cluster", o2HistType::kTH1F, {{400, 0, 100}}); + mHistManager.add("clusterE_SimpleBinning", "Energy of cluster", o2HistType::kTH1F, {{2000, 0, 200}}); mHistManager.add("clusterEtaPhi", "Eta and phi of cluster", o2HistType::kTH2F, {{100, -1, 1}, {100, 0, 2 * TMath::Pi()}}); mHistManager.add("clusterM02", "M02 of cluster", o2HistType::kTH1F, {{400, 0, 5}}); mHistManager.add("clusterM20", "M20 of cluster", o2HistType::kTH1F, {{400, 0, 2.5}}); @@ -135,6 +130,12 @@ struct ClusterMonitor { } } } + + // define cluster filter. It selects only those clusters which are of the type + // sadly passing of the string at runtime is not possible for technical region so cluster definition is + // an integer instead + Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == mClusterDefinition); + /// \brief Process EMCAL clusters that are matched to a collisions void processCollisions(collisionEvSelIt const& theCollision, selectedClusters const& clusters, o2::aod::BCs const& bcs) { @@ -144,12 +145,7 @@ struct ClusterMonitor { // currently the event selection is hard coded to kINT7 // but other selections are possible that are defined in TriggerAliases.h if (mDoEventSel && (!theCollision.alias()[kINT7])) { - LOG(debug) << "Event not selected becaus it is not kINT7, skipping"; - return; - } - mHistManager.fill(HIST("eventVertexZAll"), theCollision.posZ()); - if (mVertexCut > 0 && TMath::Abs(theCollision.posZ()) > mVertexCut) { - LOG(debug) << "Event not selected because of z-vertex cut z= " << theCollision.posZ() << " > " << mVertexCut << " cm, skipping"; + LOG(debug) << "Event not selected because it is not kINT7, skipping"; return; } mHistManager.fill(HIST("eventVertexZAll"), theCollision.posZ()); From 0b365fef771eae30765e1f944c77ab0ac9770c84 Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Fri, 29 Jul 2022 11:42:49 +0200 Subject: [PATCH 2/4] [EMCAL-689] added EMC cluster track matching and link of clusters to cells --- PWGJE/Core/JetUtilities.h | 105 ++++++++++++++++++++ PWGJE/DataModel/EMCALClusters.h | 23 ++++- PWGJE/TableProducer/emcalCorrectionTask.cxx | 73 +++++++++++++- PWGJE/Tasks/emcclustermonitor.cxx | 56 ++++++++++- 4 files changed, 253 insertions(+), 4 deletions(-) diff --git a/PWGJE/Core/JetUtilities.h b/PWGJE/Core/JetUtilities.h index 12d365257ae..212db69c30b 100644 --- a/PWGJE/Core/JetUtilities.h +++ b/PWGJE/Core/JetUtilities.h @@ -278,6 +278,111 @@ std::tuple, std::vector> MatchJetsGeometrically( return std::make_tuple(baseToTagMap, tagToBaseMap); } +/** + * Match clusters and tracks. + * + * Match cluster with tracks, where maxNumberMatches are considered in dR=maxMatchingDistance. + * If no unique match was found for a jet, an index of -1 is stored. + * The same map is created for clusters matched to tracks e.g. for electron analyses. + * + * @param clusterPhi cluster collection phi. + * @param clusterEta cluster collection eta. + * @param trackPhi track collection phi. + * @param trackEta track collection eta. + * @param maxMatchingDistance Maximum matching distance. + * @param maxNumberMatches Maximum number of matches (e.g. 5 closest). + * + * @returns (cluster to track index map, track to cluster index map) + */ +template +std::tuple>, std::vector>> MatchClustersAndTracks( + std::vector& clusterPhi, + std::vector& clusterEta, + std::vector& trackPhi, + std::vector& trackEta, + double maxMatchingDistance, + int maxNumberMatches) +{ + // test + // Validation + const std::size_t nClusters = clusterEta.size(); + const std::size_t nTracks = trackEta.size(); + if (!(nClusters && nTracks)) { + // There are no jets, so nothing to be done. + return std::make_tuple(std::vector>(nClusters, std::vector(maxNumberMatches, -1)), std::vector>(nTracks, std::vector(maxNumberMatches, -1))); + } + // Input sizes must match + if (clusterPhi.size() != clusterEta.size()) { + throw std::invalid_argument("cluster collection eta and phi sizes don't match. Check the inputs."); + } + if (trackPhi.size() != trackEta.size()) { + throw std::invalid_argument("track collection eta and phi sizes don't match. Check the inputs."); + } + + // for (std::size_t iTrack = 0; iTrack < nTracks; iTrack++) { + // if (trackEta[iTrack] == 0) + // LOG(warning) << "Track eta is 0!"; + // } + + // Build the KD-trees using vectors + // We build two trees: + // treeBase, which contains the base collection. + // treeTag, which contains the tag collection. + // The trees are built to match in two dimensions (eta, phi) + TKDTree treeCluster(clusterEta.size(), 2, 1), treeTrack(trackEta.size(), 2, 1); + // By utilizing SetData, we can avoid having to copy the data again. + treeCluster.SetData(0, clusterEta.data()); + treeCluster.SetData(1, clusterPhi.data()); + treeCluster.Build(); + treeTrack.SetData(0, trackEta.data()); + treeTrack.SetData(1, trackPhi.data()); + treeTrack.Build(); + + // Storage for the cluster matching indices. + std::vector> matchIndexTrack(nClusters, std::vector(maxNumberMatches, -1)); + std::vector> matchIndexCluster(nTracks, std::vector(maxNumberMatches, -1)); + + // Find the track closest to each cluster. + for (std::size_t iCluster = 0; iCluster < nClusters; iCluster++) { + T point[2] = {clusterEta[iCluster], clusterPhi[iCluster]}; + int index[50]; // size 50 for safety + T distance[50]; // size 50 for safery + std::fill_n(index, 50, -1); + std::fill_n(distance, 50, std::numeric_limits::max()); + treeTrack.FindNearestNeighbors(point, maxNumberMatches, index, distance); + // test whether indices are matching: + matchIndexTrack[iCluster] = std::vector(maxNumberMatches); + for (int m = 0; m < maxNumberMatches; m++) { + if (index[m] >= 0 && distance[m] < maxMatchingDistance) { + matchIndexTrack[iCluster][m] = index[m]; + } else { + // no match or no more matches found, fill -1 + matchIndexTrack[iCluster][m] = -1; + } + } + } + + // Find the base jet closest to each tag jet + for (std::size_t iTrack = 0; iTrack < nTracks; iTrack++) { + T point[2] = {trackEta[iTrack], trackPhi[iTrack]}; + int index[50]; // size 50 for safety + T distance[50]; // size 50 for safery + std::fill_n(index, 50, -1); + std::fill_n(distance, 50, std::numeric_limits::max()); + treeCluster.FindNearestNeighbors(point, maxNumberMatches, index, distance); + matchIndexCluster[iTrack] = std::vector(maxNumberMatches); + // loop over maxNumberMatches closest matches + for (int m = 0; m < maxNumberMatches; m++) { + if (index[m] >= 0 && distance[m] < maxMatchingDistance) { + matchIndexCluster[iTrack][m] = index[m]; + } else { + // no match jet or no more matches found, fill -1 + matchIndexCluster[iTrack][m] = -1; + } + } + } + return std::make_tuple(matchIndexTrack, matchIndexCluster); +} }; // namespace JetUtilities #endif diff --git a/PWGJE/DataModel/EMCALClusters.h b/PWGJE/DataModel/EMCALClusters.h index 8dc8f5d18c7..d62222dbab0 100644 --- a/PWGJE/DataModel/EMCALClusters.h +++ b/PWGJE/DataModel/EMCALClusters.h @@ -89,6 +89,27 @@ DECLARE_SOA_TABLE(EMCALAmbiguousClusters, "AOD", "EMCALAMBCLUS", //! using EMCALCluster = EMCALClusters::iterator; using EMCALAmbiguousCluster = EMCALAmbiguousClusters::iterator; -} // namespace o2::aod +namespace emcalclustercell +{ +// declare index column pointing to cluster table +DECLARE_SOA_INDEX_COLUMN(EMCALCluster, emcalcluster); //! linked to EMCalClusters table +DECLARE_SOA_INDEX_COLUMN(Calo, calo); //! linked to calo cells +// declare index column pointing to ambiguous cluster table +DECLARE_SOA_INDEX_COLUMN(EMCALAmbiguousCluster, emcalambiguouscluster); //! linked to EMCalAmbiguousClusters table +} // namespace emcalclustercell +DECLARE_SOA_TABLE(EMCALClusterCells, "AOD", "EMCCLUSCELLS", //! + o2::soa::Index<>, emcalclustercell::EMCALClusterId, emcalclustercell::CaloId); //! +DECLARE_SOA_TABLE(EMCALAmbiguousClusterCells, "AOD", "EMCAMBBCLUSCLS", //! + o2::soa::Index<>, emcalclustercell::EMCALAmbiguousClusterId, emcalclustercell::CaloId); //! +using EMCALClusterCell = EMCALClusterCells::iterator; +using EMCALAmbiguousClusterCell = EMCALAmbiguousClusterCells::iterator; +namespace emcalmatchedtrack +{ +DECLARE_SOA_INDEX_COLUMN(Track, track); //! linked to Track table only for tracks that were matched +} // namespace emcalmatchedtrack +DECLARE_SOA_TABLE(EMCALMatchedTracks, "AOD", "EMCMATCHTRACKS", //! + o2::soa::Index<>, emcalclustercell::EMCALClusterId, emcalmatchedtrack::TrackId); //! +using EMCALMatchedTrack = EMCALMatchedTracks::iterator; +} // namespace o2::aod #endif diff --git a/PWGJE/TableProducer/emcalCorrectionTask.cxx b/PWGJE/TableProducer/emcalCorrectionTask.cxx index f82060f53d9..2d004e48417 100644 --- a/PWGJE/TableProducer/emcalCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalCorrectionTask.cxx @@ -32,6 +32,7 @@ #include "EMCALBase/Geometry.h" #include "EMCALBase/ClusterFactory.h" #include "EMCALReconstruction/Clusterizer.h" +#include "PWGJE/Core/JetUtilities.h" #include "TVector2.h" using namespace o2; @@ -41,11 +42,17 @@ using namespace o2::framework::expressions; struct EmcalCorrectionTask { Produces clusters; Produces clustersAmbiguous; + Produces clustercells; // cells belonging to given cluster + Produces clustercellsambiguous; + Produces matchedTracks; + Preslice perCollision = aod::track::collisionId; // Options for the clusterization // 1 corresponds to EMCAL cells based on the Run2 definition. Configurable selectedCellType{"selectedCellType", 1, "EMCAL Cell type"}; Configurable clusterDefinitions{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default. Multiple definitions can be specified separated by comma"}; + Configurable maxMatchingDistance{"maxMatchingDistance", 0.4f, "Max matching distance track-cluster"}; + // CDB service (for geometry) Service mCcdbManager; @@ -56,6 +63,8 @@ struct EmcalCorrectionTask { std::vector>> mClusterFactories; // Cells and clusters std::vector mEmcalCells; + // map of cellId (local in BC) to global cell index in cell table in AO2D + std::map mCellIdToCellGlobalIndex; std::vector mAnalysisClusters; std::vector mClusterDefinitions; @@ -135,15 +144,17 @@ struct EmcalCorrectionTask { // void process(aod::BCs const& bcs, aod::Collision const& collision, aod::Calos const& cells) // Appears to need the BC to be accessed to be available in the collision table... - void process(aod::BC const& bc, aod::Collisions const& collisions, aod::Calos const& cells) + void process(aod::BC const& bc, aod::Collisions const& collisions, aod::Tracks const& tracks, aod::Calos const& cells) { LOG(debug) << "Starting process."; // Convert aod::Calo to o2::emcal::Cell which can be used with the clusterizer. // In particular, we need to filter only EMCAL cells. mEmcalCells.clear(); + mCellIdToCellGlobalIndex.clear(); + int c = 0; for (auto& cell : cells) { if (cell.caloType() != selectedCellType) { - // LOG(debug) << "Rejected"; + LOG(debug) << "Rejected"; continue; } // LOG(debug) << "Cell E: " << cell.getEnergy(); @@ -153,6 +164,9 @@ struct EmcalCorrectionTask { cell.amplitude(), cell.time(), o2::emcal::intToChannelType(cell.cellType()))); + mCellIdToCellGlobalIndex.insert(std::make_pair(c, cell.globalIndex())); + LOG(debug) << "Creating map " << c << " -> " << cell.globalIndex(); + c++; } LOG(debug) << "Number of cells (CF): " << mEmcalCells.size(); @@ -215,8 +229,41 @@ struct EmcalCorrectionTask { vz = col.posZ(); hasCollision = true; + // store positions of all tracks of collision + auto groupedTracks = tracks.sliceBy(perCollision, col.globalIndex()); + + std::vector trackPhi; + std::vector trackEta; + std::vector trackGlobalIndex; + for (auto& track : groupedTracks) { + // TODO this actually needs to use the eta phi + // of track propagated to EMC surface! Will be provided centrally according to Ruben + // TODO only consider tracks in current emcal/dcal acceptanc + trackPhi.emplace_back(TVector2::Phi_0_2pi(track.phi())); + trackEta.emplace_back(track.eta()); + trackGlobalIndex.emplace_back(track.globalIndex()); + } + + std::vector clusterPhi; + std::vector clusterEta; + + // TODO one loop that could in principle be combined with the other loop to improve performance + for (const auto& cluster : mAnalysisClusters) { + // Determine the cluster eta, phi, correcting for the vertex position. + auto pos = cluster.getGlobalPosition(); + pos = pos - math_utils::Point3D{vx, vy, vz}; + // Normalize the vector and rescale by energy. + pos *= (cluster.E() / std::sqrt(pos.Mag2())); + clusterPhi.emplace_back(TVector2::Phi_0_2pi(pos.Phi())); + clusterEta.emplace_back(pos.Eta()); + } + auto&& [clusterToTrackIndexMap, trackToClusterIndexMap] = JetUtilities::MatchClustersAndTracks(clusterPhi, clusterEta, trackPhi, trackEta, maxMatchingDistance, 5); // we found a collision, put the clusters into the none ambiguous table clusters.reserve(mAnalysisClusters.size()); + int cellindex = -1; + float cellampfraction = -1; + + unsigned int k = 0; for (const auto& cluster : mAnalysisClusters) { // Determine the cluster eta, phi, correcting for the vertex position. @@ -231,9 +278,23 @@ struct EmcalCorrectionTask { cluster.getM02(), cluster.getM20(), cluster.getNCells(), cluster.getClusterTime(), cluster.getIsExotic(), cluster.getDistanceToBadChannel(), cluster.getNExMax(), static_cast(mClusterDefinitions.at(i))); + clustercells.reserve(cluster.getNCells()); + // loop over cells in cluster and save to table + for (int ncell = 0; ncell < cluster.getNCells(); ncell++) { + cellindex = cluster.getCellIndex(ncell); + LOG(debug) << "trying to find cell index " << cellindex << " in map"; + clustercells(clusters.lastIndex(), mCellIdToCellGlobalIndex.at(cellindex)); + } // fill histograms hClusterE->Fill(cluster.E()); hClusterEtaPhi->Fill(pos.Eta(), TVector2::Phi_0_2pi(pos.Phi())); + for (int iTrack = 0; iTrack < clusterToTrackIndexMap[k].size(); iTrack++) { + if (clusterToTrackIndexMap[k][iTrack] >= 0) { + LOG(debug) << "Found track " << trackGlobalIndex[clusterToTrackIndexMap[k][iTrack]] << " in cluster " << cluster.getID(); + matchedTracks(clusters.lastIndex(), trackGlobalIndex[clusterToTrackIndexMap[k][iTrack]]); + } + } + k++; } // end of cluster loop } // end of collision loop } @@ -244,6 +305,8 @@ struct EmcalCorrectionTask { // Store the clusters in the table where a mathcing collision could // be identified. if (!hasCollision) { // ambiguous + int cellindex = -1; + float cellampfraction = -1; clustersAmbiguous.reserve(mAnalysisClusters.size()); for (const auto& cluster : mAnalysisClusters) { auto pos = cluster.getGlobalPosition(); @@ -257,6 +320,12 @@ struct EmcalCorrectionTask { clustersAmbiguous(bc, cluster.getID(), cluster.E(), cluster.getCoreEnergy(), pos.Eta(), TVector2::Phi_0_2pi(pos.Phi()), cluster.getM02(), cluster.getM20(), cluster.getNCells(), cluster.getClusterTime(), cluster.getIsExotic(), cluster.getDistanceToBadChannel(), cluster.getNExMax(), static_cast(mClusterDefinitions.at(i))); + clustercellsambiguous.reserve(cluster.getNCells()); + for (int ncell = 0; ncell < cluster.getNCells(); ncell++) { + cellindex = cluster.getCellIndex(ncell); + cellampfraction = cluster.getCellAmplitudeFraction(ncell); + clustercellsambiguous(clustersAmbiguous.lastIndex(), mCellIdToCellGlobalIndex.at(cellindex)); + } } } LOG(debug) << "Cluster loop done for clusterizer " << i; diff --git a/PWGJE/Tasks/emcclustermonitor.cxx b/PWGJE/Tasks/emcclustermonitor.cxx index 03d3f7e5f6f..ae11266c816 100644 --- a/PWGJE/Tasks/emcclustermonitor.cxx +++ b/PWGJE/Tasks/emcclustermonitor.cxx @@ -59,6 +59,9 @@ struct ClusterMonitor { HistogramRegistry mHistManager{"ClusterMonitorHistograms"}; o2::emcal::Geometry* mGeometry = nullptr; + Preslice perCluster = o2::aod::emcalclustercell::emcalclusterId; + Preslice perClusterAmb = o2::aod::emcalclustercell::emcalclusterId; + Preslice perClusterMatchedTracks = o2::aod::emcalclustercell::emcalclusterId; // configurable parameters // TODO adapt mDoEventSel switch to also allow selection of other triggers (e.g. EMC7) Configurable mDoEventSel{"doEventSel", 0, "demand kINT7"}; @@ -100,6 +103,7 @@ struct ClusterMonitor { // cluster properties (matched clusters) mHistManager.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis}); + mHistManager.add("clusterEMatched", "Energy of cluster (with match)", o2HistType::kTH1F, {energyAxis}); mHistManager.add("clusterE_SimpleBinning", "Energy of cluster", o2HistType::kTH1F, {{2000, 0, 200}}); mHistManager.add("clusterEtaPhi", "Eta and phi of cluster", o2HistType::kTH2F, {{100, -1, 1}, {100, 0, 2 * TMath::Pi()}}); mHistManager.add("clusterM02", "M02 of cluster", o2HistType::kTH1F, {{400, 0, 5}}); @@ -108,6 +112,9 @@ struct ClusterMonitor { mHistManager.add("clusterNCells", "Number of cells in cluster", o2HistType::kTH1I, {{50, 0, 50}}); mHistManager.add("clusterDistanceToBadChannel", "Distance to bad channel", o2HistType::kTH1F, {{100, 0, 100}}); mHistManager.add("clusterTimeVsE", "Cluster time vs energy", o2HistType::kTH2F, {timeAxis, energyAxis}); + mHistManager.add("clusterAmpFractionLeadingCell", "Fraction of energy in leading cell", o2HistType::kTH1F, {{100, 0, 1}}); + mHistManager.add("clusterTM_dEtadPhi", "cluster trackmatching dEta/dPhi", o2HistType::kTH2F, {{100, -0.4, 0.4}, {100, -0.4, 0.4}}); + mHistManager.add("clusterTM_EoverP_E", "cluster E/p (dEtadPhi<0.05)", o2HistType::kTH2F, {{500, 0, 10}, {200, 0, 100}}); if (mVetoBCID->length()) { std::stringstream parser(mVetoBCID.value); @@ -137,7 +144,7 @@ struct ClusterMonitor { Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == mClusterDefinition); /// \brief Process EMCAL clusters that are matched to a collisions - void processCollisions(collisionEvSelIt const& theCollision, selectedClusters const& clusters, o2::aod::BCs const& bcs) + void processCollisions(collisionEvSelIt const& theCollision, selectedClusters const& clusters, o2::aod::EMCALClusterCells const& emccluscells, o2::aod::Calos const& allcalos, o2::aod::EMCALMatchedTracks const& matchedtracks, o2::aod::Tracks const& alltrack) { mHistManager.fill(HIST("eventsAll"), 1); @@ -156,6 +163,7 @@ struct ClusterMonitor { mHistManager.fill(HIST("eventsSelected"), 1); mHistManager.fill(HIST("eventVertexZSelected"), theCollision.posZ()); + LOG(debug) << "bunch crossing ID" << theCollision.bcId(); // loop over all clusters from accepted collision // auto eventClusters = clusters.select(o2::aod::emcalcluster::bcId == theCollision.bc().globalBC()); for (const auto& cluster : clusters) { @@ -176,6 +184,52 @@ struct ClusterMonitor { mHistManager.fill(HIST("clusterNLM"), cluster.nlm()); mHistManager.fill(HIST("clusterNCells"), cluster.nCells()); mHistManager.fill(HIST("clusterDistanceToBadChannel"), cluster.distanceToBadChannel()); + // loop over cells in cluster + LOG(debug) << "Cluster energy: " << cluster.energy(); + LOG(debug) << "Cluster index: " << cluster.index(); + LOG(debug) << "ncells in cluster: " << cluster.nCells(); + LOG(debug) << "real cluster id" << cluster.id(); + + // example of loop over all cells of current cluster + // cell.calo() allows to access the cell properties as defined in AnalysisDataModel + // In this exammple, we loop over all cells and find the cell of maximum energy and plot the fraction + // it carries of the whole cluster + auto cellsofcluster = emccluscells.sliceBy(perCluster, cluster.globalIndex()); + double maxamp = 0; + double ampfraction = 0; + for (const auto& cell : cellsofcluster) { + // example how to get any information of the cell associated with cluster + LOG(debug) << "Cell ID:" << cell.calo().amplitude() << " Time " << cell.calo().time(); + if (cell.calo().amplitude() > maxamp) { + maxamp = cell.calo().amplitude(); + } + } + ampfraction = maxamp / cluster.energy(); + mHistManager.fill(HIST("clusterAmpFractionLeadingCell"), ampfraction); + + // Example of loop over tracks matched to cluster within dR=0.4, where only the + // 5 most closest tracks are stored. If needed the number of tracks can be later + // increased in the correction framework. Access to all track properties via match.track() + // In this example the counter t is just used to only look at the closest match + double dEta, dPhi, dR; + auto tracksofcluster = matchedtracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex()); + int t = 0; + for (const auto& match : tracksofcluster) { + // exmple of how to access any property of the matched tracks (tracks are sorted by how close they are to cluster) + LOG(debug) << "Pt of match" << match.track().pt(); + // only consider closest match + if (t == 0) { + dEta = match.track().eta() - cluster.eta(); + dPhi = match.track().phi() - cluster.phi(); + dR = TMath::Sqrt(dEta * dEta + dPhi * dPhi); + if (dEta < 0.05 && dPhi < 0.05) { + mHistManager.fill(HIST("clusterTM_EoverP_E"), cluster.energy() / match.track().p(), cluster.energy()); + } + mHistManager.fill(HIST("clusterTM_dEtadPhi"), dEta, dPhi); + } + + t++; + } } } PROCESS_SWITCH(ClusterMonitor, processCollisions, "Process clusters from collision", false); From e9fc3a62d33ff17172b8352277dd2b08caccf2cb Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Fri, 29 Jul 2022 12:06:07 +0200 Subject: [PATCH 3/4] [EMCAL-689] fixed trailing whitespace --- PWGJE/Core/JetUtilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGJE/Core/JetUtilities.h b/PWGJE/Core/JetUtilities.h index 212db69c30b..2f3d8c2adce 100644 --- a/PWGJE/Core/JetUtilities.h +++ b/PWGJE/Core/JetUtilities.h @@ -281,7 +281,7 @@ std::tuple, std::vector> MatchJetsGeometrically( /** * Match clusters and tracks. * - * Match cluster with tracks, where maxNumberMatches are considered in dR=maxMatchingDistance. + * Match cluster with tracks, where maxNumberMatches are considered in dR=maxMatchingDistance. * If no unique match was found for a jet, an index of -1 is stored. * The same map is created for clusters matched to tracks e.g. for electron analyses. * From a089ca488f2660fc793c0677ee5e8155cd381544 Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Fri, 29 Jul 2022 12:26:18 +0200 Subject: [PATCH 4/4] [EMCAL-689] fixed tcomiler warnings --- PWGJE/TableProducer/emcalCorrectionTask.cxx | 5 +---- PWGJE/Tasks/emcclustermonitor.cxx | 3 +-- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/PWGJE/TableProducer/emcalCorrectionTask.cxx b/PWGJE/TableProducer/emcalCorrectionTask.cxx index 2d004e48417..2189910fd1e 100644 --- a/PWGJE/TableProducer/emcalCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalCorrectionTask.cxx @@ -261,7 +261,6 @@ struct EmcalCorrectionTask { // we found a collision, put the clusters into the none ambiguous table clusters.reserve(mAnalysisClusters.size()); int cellindex = -1; - float cellampfraction = -1; unsigned int k = 0; for (const auto& cluster : mAnalysisClusters) { @@ -288,7 +287,7 @@ struct EmcalCorrectionTask { // fill histograms hClusterE->Fill(cluster.E()); hClusterEtaPhi->Fill(pos.Eta(), TVector2::Phi_0_2pi(pos.Phi())); - for (int iTrack = 0; iTrack < clusterToTrackIndexMap[k].size(); iTrack++) { + for (unsigned int iTrack = 0; iTrack < clusterToTrackIndexMap[k].size(); iTrack++) { if (clusterToTrackIndexMap[k][iTrack] >= 0) { LOG(debug) << "Found track " << trackGlobalIndex[clusterToTrackIndexMap[k][iTrack]] << " in cluster " << cluster.getID(); matchedTracks(clusters.lastIndex(), trackGlobalIndex[clusterToTrackIndexMap[k][iTrack]]); @@ -306,7 +305,6 @@ struct EmcalCorrectionTask { // be identified. if (!hasCollision) { // ambiguous int cellindex = -1; - float cellampfraction = -1; clustersAmbiguous.reserve(mAnalysisClusters.size()); for (const auto& cluster : mAnalysisClusters) { auto pos = cluster.getGlobalPosition(); @@ -323,7 +321,6 @@ struct EmcalCorrectionTask { clustercellsambiguous.reserve(cluster.getNCells()); for (int ncell = 0; ncell < cluster.getNCells(); ncell++) { cellindex = cluster.getCellIndex(ncell); - cellampfraction = cluster.getCellAmplitudeFraction(ncell); clustercellsambiguous(clustersAmbiguous.lastIndex(), mCellIdToCellGlobalIndex.at(cellindex)); } } diff --git a/PWGJE/Tasks/emcclustermonitor.cxx b/PWGJE/Tasks/emcclustermonitor.cxx index ae11266c816..6b7d49a8e2c 100644 --- a/PWGJE/Tasks/emcclustermonitor.cxx +++ b/PWGJE/Tasks/emcclustermonitor.cxx @@ -211,7 +211,7 @@ struct ClusterMonitor { // 5 most closest tracks are stored. If needed the number of tracks can be later // increased in the correction framework. Access to all track properties via match.track() // In this example the counter t is just used to only look at the closest match - double dEta, dPhi, dR; + double dEta, dPhi; auto tracksofcluster = matchedtracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex()); int t = 0; for (const auto& match : tracksofcluster) { @@ -221,7 +221,6 @@ struct ClusterMonitor { if (t == 0) { dEta = match.track().eta() - cluster.eta(); dPhi = match.track().phi() - cluster.phi(); - dR = TMath::Sqrt(dEta * dEta + dPhi * dPhi); if (dEta < 0.05 && dPhi < 0.05) { mHistManager.fill(HIST("clusterTM_EoverP_E"), cluster.energy() / match.track().p(), cluster.energy()); }