diff --git a/Analysis/DataModel/include/AnalysisDataModel/EMCALClusters.h b/Analysis/DataModel/include/AnalysisDataModel/EMCALClusters.h new file mode 100644 index 0000000000000..a75ad1f8ca00b --- /dev/null +++ b/Analysis/DataModel/include/AnalysisDataModel/EMCALClusters.h @@ -0,0 +1,36 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// Table definitions for EMCAL analysis clusters +// +// Author: Raymond Ehlers + +#pragma once + +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +namespace emcalcluster +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); +DECLARE_SOA_COLUMN(Energy, energy, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(M02, m02, float); +} // namespace emcalcluster + +DECLARE_SOA_TABLE(EMCALClusters, "AOD", "EMCALCLUSTERS", + o2::soa::Index<>, emcalcluster::CollisionId, emcalcluster::Energy, + emcalcluster::Eta, emcalcluster::Phi, emcalcluster::M02); + +using EMCALCluster = EMCALClusters::iterator; + +} // namespace o2::aod diff --git a/Analysis/Tasks/CMakeLists.txt b/Analysis/Tasks/CMakeLists.txt index 662a3090d1ec0..3381c8457c3ef 100644 --- a/Analysis/Tasks/CMakeLists.txt +++ b/Analysis/Tasks/CMakeLists.txt @@ -87,3 +87,8 @@ o2_add_dpl_workflow(weak-decay-indices SOURCES weakDecayIndices.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(emcal-correction-task + SOURCES emcalCorrectionTask.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsBase O2::EMCALBase O2::EMCALReconstruction + COMPONENT_NAME Analysis) diff --git a/Analysis/Tasks/emcalCorrectionTask.cxx b/Analysis/Tasks/emcalCorrectionTask.cxx new file mode 100644 index 0000000000000..dfba6072ae069 --- /dev/null +++ b/Analysis/Tasks/emcalCorrectionTask.cxx @@ -0,0 +1,198 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// EMCAL Correction Task +// +// Author: Raymond Ehlers + +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" + +#include "DetectorsBase/GeometryManager.h" + +#include "AnalysisDataModel/EMCALClusters.h" + +#include "DataFormatsEMCAL/Cell.h" +#include "DataFormatsEMCAL/Constants.h" +#include "DataFormatsEMCAL/AnalysisCluster.h" +#include "EMCALBase/Geometry.h" +#include "EMCALBase/ClusterFactory.h" +#include "EMCALReconstruction/Clusterizer.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct EmcalCorrectionTask { + Produces clusters; + + // Options for the clusterization + // 1 corresponds to EMCAL cells based on the Run2 definition. + Configurable selectedCellType{"selectedCellType", 1, "EMCAL Cell type"}; + Configurable seedEnergy{"seedEnergy", 0.1, "Clusterizer seed energy."}; + Configurable minCellEnergy{"minCellEnergy", 0.05, "Clusterizer minimum cell energy."}; + // TODO: Check this range, especially after change to the conversion... + Configurable timeCut{"timeCut", 10000, "Cell time cut"}; + Configurable timeMin{"timeMin", 0, "Min cell time"}; + Configurable timeMax{"timeMax", 10000, "Max cell time"}; + Configurable enableEnergyGradientCut{"enableEnergyGradientCut", true, "Enable energy gradient cut."}; + Configurable gradientCut{"gradientCut", 0.03, "Clusterizer energy gradient cut."}; + + // Clusterizer and related + // Apparently streaming these objects really doesn't work, and causes problems for setting up the workflow. + // So we use unique_ptr and define them below. + std::unique_ptr> mClusterizer; + std::unique_ptr> mClusterFactory; + // Cells and clusters + std::vector mEmcalCells; + std::vector mAnalysisClusters; + + // QA + // NOTE: This is not comprehensive. + OutputObj hCellE{"hCellE"}; + OutputObj hCellTowerID{"hCellTowerID"}; + OutputObj hCellEtaPhi{"hCellEtaPhi"}; + OutputObj hCellRowCol{"hCellRowCol"}; + OutputObj hClusterE{"hClusterE"}; + OutputObj hClusterEtaPhi{"hClusterEtaPhi"}; + + void init(InitContext const&) + { + LOG(DEBUG) << "Start init!"; + // NOTE: The geometry manager isn't necessary just to load the EMCAL geometry. + // However, it _is_ necessary for loading the misalignment matrices as of September 2020 + // Eventually, those matrices will be moved to the CCDB, but it's not yet ready. + // FIXME: Hardcoded for run 2 + o2::base::GeometryManager::loadGeometry(); // for generating full clusters + LOG(DEBUG) << "After load geometry!"; + o2::emcal::Geometry* geometry = o2::emcal::Geometry::GetInstanceFromRunNumber(223409); + if (!geometry) { + LOG(ERROR) << "Failure accessing geometry"; + } + + // Setup clusterizer + LOG(DEBUG) << "Init clusterizer!"; + mClusterizer = decltype(mClusterizer)(new o2::emcal::Clusterizer()); + mClusterizer->initialize(timeCut, timeMin, timeMax, gradientCut, enableEnergyGradientCut, seedEnergy, minCellEnergy); + mClusterizer->setGeometry(geometry); + LOG(DEBUG) << "Done with clusterizer. Setup cluster factory."; + // Setup cluster factory. + mClusterFactory = decltype(mClusterFactory)(new o2::emcal::ClusterFactory()); + LOG(DEBUG) << "Completed init!"; + + // Setup QA hists. + hCellE.setObject(new TH1F("hCellE", "hCellE", 200, 0.0, 100)); + hCellTowerID.setObject(new TH1I("hCellTowerID", "hCellTowerID", 20000, 0, 20000)); + hCellEtaPhi.setObject(new TH2F("hCellEtaPhi", "hCellEtaPhi", 160, -0.8, 0.8, 72, 0, 2 * 3.14159)); + // NOTE: Reversed column and row because it's more natural for presentatin. + hCellRowCol.setObject(new TH2I("hCellRowCol", "hCellRowCol;Column;Row", 97, 0, 97, 600, 0, 600)); + hClusterE.setObject(new TH1F("hClusterE", "hClusterE", 200, 0.0, 100)); + hClusterEtaPhi.setObject(new TH2F("hClusterEtaPhi", "hClusterEtaPhi", 160, -0.8, 0.8, 72, 0, 2 * 3.14159)); + } + + //void process(aod::Collision const& collision, soa::Filtered const& fullTracks, aod::Calos const& cells) + //void process(aod::Collision const& collision, aod::Tracks const& tracks, aod::Calos const& cells) + //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::Collision const& collision, aod::Calos const& cells, aod::BCs const& bcs) + { + 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(); + for (auto& cell : cells) { + if (cell.caloType() != selectedCellType || cell.bc() != collision.bc()) { + //LOG(DEBUG) << "Rejected"; + continue; + } + //LOG(DEBUG) << "Cell E: " << cell.getEnergy(); + //LOG(DEBUG) << "Cell E: " << cell; + + mEmcalCells.emplace_back(o2::emcal::Cell( + cell.cellNumber(), + cell.amplitude(), + cell.time(), + o2::emcal::intToChannelType(cell.cellType()))); + } + + // Cell QA + // For convenience, use the clusterizer stored geometry to get the eta-phi + for (auto& cell : mEmcalCells) { + hCellE->Fill(cell.getEnergy()); + hCellTowerID->Fill(cell.getTower()); + auto res = mClusterizer->getGeometry()->EtaPhiFromIndex(cell.getTower()); + hCellEtaPhi->Fill(std::get<0>(res), std::get<1>(res)); + res = mClusterizer->getGeometry()->GlobalRowColFromIndex(cell.getTower()); + // NOTE: Reversed column and row because it's more natural for presentatin. + hCellRowCol->Fill(std::get<1>(res), std::get<0>(res)); + } + + // TODO: Helpful for now, but should be removed. + LOG(DEBUG) << "Converted EMCAL cells"; + for (auto& cell : mEmcalCells) { + LOG(DEBUG) << cell.getTower() << ": E: " << cell.getEnergy() << ", time: " << cell.getTimeStamp() << ", type: " << cell.getType(); + } + + LOG(INFO) << "Converted cells. Contains: " << mEmcalCells.size() << ". Originally " << cells.size() << ". About to run clusterizer."; + + // Run the clusterizer + mClusterizer->findClusters(mEmcalCells); + LOG(DEBUG) << "Found clusters."; + auto emcalClusters = mClusterizer->getFoundClusters(); + auto emcalClustersInputIndices = mClusterizer->getFoundClustersInputIndices(); + LOG(DEBUG) << "Retrieved results. About to setup cluster factory."; + + // Convert to analysis clusters. + // First, the cluster factory requires cluster and cell information in order to build the clusters. + mAnalysisClusters.clear(); + mClusterFactory->reset(); + mClusterFactory->setClustersContainer(*emcalClusters); + mClusterFactory->setCellsContainer(mEmcalCells); + mClusterFactory->setCellsIndicesContainer(*emcalClustersInputIndices); + LOG(DEBUG) << "Cluster factory set up."; + + // Convert to analysis clusters. + for (int icl = 0; icl < mClusterFactory->getNumberOfClusters(); icl++) { + auto analysisCluster = mClusterFactory->buildCluster(icl); + mAnalysisClusters.emplace_back(analysisCluster); + } + LOG(DEBUG) << "Converted to analysis clusters."; + + // Store the clusters in the table + clusters.reserve(mAnalysisClusters.size()); + for (const auto& cluster : mAnalysisClusters) { + // Determine the cluster eta, phi, correcting for the vertex position. + auto pos = cluster.getGlobalPosition(); + pos = pos - math_utils::Point3D{collision.posX(), collision.posY(), collision.posZ()}; + // Normalize the vector and rescale by energy. + pos /= (cluster.E() / std::sqrt(pos.Mag2())); + + // We have our necessary properties. Now we store outputs + //LOG(DEBUG) << "Cluster E: " << cluster.E(); + clusters(collision, cluster.E(), pos.Eta(), pos.Phi(), cluster.getM02()); + //if (cluster.E() < 0.300) { + // continue; + //} + hClusterE->Fill(cluster.E()); + hClusterEtaPhi->Fill(pos.Eta(), pos.Phi()); + } + LOG(DEBUG) << "Done with process."; + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("emcal-correction-task")}; +}