Skip to content

Commit 7ed4626

Browse files
raymondEhlersraymondEhlers
authored
Initial version of EMCAL Correction Task (#5074)
1 parent 585bc0f commit 7ed4626

File tree

3 files changed

+239
-0
lines changed

3 files changed

+239
-0
lines changed

‎Analysis/DataModel/include/AnalysisDataModel/EMCALClusters.h‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
// Table definitions for EMCAL analysis clusters
12+
//
13+
// Author: Raymond Ehlers
14+
15+
#pragma once
16+
17+
#include "Framework/AnalysisDataModel.h"
18+
19+
namespace o2::aod
20+
{
21+
namespace emcalcluster
22+
{
23+
DECLARE_SOA_INDEX_COLUMN(Collision, collision);
24+
DECLARE_SOA_COLUMN(Energy, energy, float);
25+
DECLARE_SOA_COLUMN(Eta, eta, float);
26+
DECLARE_SOA_COLUMN(Phi, phi, float);
27+
DECLARE_SOA_COLUMN(M02, m02, float);
28+
} // namespace emcalcluster
29+
30+
DECLARE_SOA_TABLE(EMCALClusters, "AOD", "EMCALCLUSTERS",
31+
o2::soa::Index<>, emcalcluster::CollisionId, emcalcluster::Energy,
32+
emcalcluster::Eta, emcalcluster::Phi, emcalcluster::M02);
33+
34+
using EMCALCluster = EMCALClusters::iterator;
35+
36+
} // namespace o2::aod

‎Analysis/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,3 +87,8 @@ o2_add_dpl_workflow(weak-decay-indices
8787
SOURCES weakDecayIndices.cxx
8888
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel
8989
COMPONENT_NAME Analysis)
90+
91+
o2_add_dpl_workflow(emcal-correction-task
92+
SOURCES emcalCorrectionTask.cxx
93+
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsBase O2::EMCALBase O2::EMCALReconstruction
94+
COMPONENT_NAME Analysis)

‎Analysis/Tasks/emcalCorrectionTask.cxx‎

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
// EMCAL Correction Task
12+
//
13+
// Author: Raymond Ehlers
14+
15+
#include <cmath>
16+
17+
#include "Framework/runDataProcessing.h"
18+
#include "Framework/AnalysisTask.h"
19+
#include "Framework/AnalysisDataModel.h"
20+
#include "Framework/ASoA.h"
21+
22+
#include "DetectorsBase/GeometryManager.h"
23+
24+
#include "AnalysisDataModel/EMCALClusters.h"
25+
26+
#include "DataFormatsEMCAL/Cell.h"
27+
#include "DataFormatsEMCAL/Constants.h"
28+
#include "DataFormatsEMCAL/AnalysisCluster.h"
29+
#include "EMCALBase/Geometry.h"
30+
#include "EMCALBase/ClusterFactory.h"
31+
#include "EMCALReconstruction/Clusterizer.h"
32+
33+
using namespace o2;
34+
using namespace o2::framework;
35+
using namespace o2::framework::expressions;
36+
37+
struct EmcalCorrectionTask {
38+
Produces<o2::aod::EMCALClusters> clusters;
39+
40+
// Options for the clusterization
41+
// 1 corresponds to EMCAL cells based on the Run2 definition.
42+
Configurable<int> selectedCellType{"selectedCellType", 1, "EMCAL Cell type"};
43+
Configurable<double> seedEnergy{"seedEnergy", 0.1, "Clusterizer seed energy."};
44+
Configurable<double> minCellEnergy{"minCellEnergy", 0.05, "Clusterizer minimum cell energy."};
45+
// TODO: Check this range, especially after change to the conversion...
46+
Configurable<double> timeCut{"timeCut", 10000, "Cell time cut"};
47+
Configurable<double> timeMin{"timeMin", 0, "Min cell time"};
48+
Configurable<double> timeMax{"timeMax", 10000, "Max cell time"};
49+
Configurable<bool> enableEnergyGradientCut{"enableEnergyGradientCut", true, "Enable energy gradient cut."};
50+
Configurable<double> gradientCut{"gradientCut", 0.03, "Clusterizer energy gradient cut."};
51+
52+
// Clusterizer and related
53+
// Apparently streaming these objects really doesn't work, and causes problems for setting up the workflow.
54+
// So we use unique_ptr and define them below.
55+
std::unique_ptr<o2::emcal::Clusterizer<o2::emcal::Cell>> mClusterizer;
56+
std::unique_ptr<o2::emcal::ClusterFactory<o2::emcal::Cell>> mClusterFactory;
57+
// Cells and clusters
58+
std::vector<o2::emcal::Cell> mEmcalCells;
59+
std::vector<o2::emcal::AnalysisCluster> mAnalysisClusters;
60+
61+
// QA
62+
// NOTE: This is not comprehensive.
63+
OutputObj<TH1F> hCellE{"hCellE"};
64+
OutputObj<TH1I> hCellTowerID{"hCellTowerID"};
65+
OutputObj<TH2F> hCellEtaPhi{"hCellEtaPhi"};
66+
OutputObj<TH2I> hCellRowCol{"hCellRowCol"};
67+
OutputObj<TH1F> hClusterE{"hClusterE"};
68+
OutputObj<TH2F> hClusterEtaPhi{"hClusterEtaPhi"};
69+
70+
void init(InitContext const&)
71+
{
72+
LOG(DEBUG) << "Start init!";
73+
// NOTE: The geometry manager isn't necessary just to load the EMCAL geometry.
74+
// However, it _is_ necessary for loading the misalignment matrices as of September 2020
75+
// Eventually, those matrices will be moved to the CCDB, but it's not yet ready.
76+
// FIXME: Hardcoded for run 2
77+
o2::base::GeometryManager::loadGeometry(); // for generating full clusters
78+
LOG(DEBUG) << "After load geometry!";
79+
o2::emcal::Geometry* geometry = o2::emcal::Geometry::GetInstanceFromRunNumber(223409);
80+
if (!geometry) {
81+
LOG(ERROR) << "Failure accessing geometry";
82+
}
83+
84+
// Setup clusterizer
85+
LOG(DEBUG) << "Init clusterizer!";
86+
mClusterizer = decltype(mClusterizer)(new o2::emcal::Clusterizer<o2::emcal::Cell>());
87+
mClusterizer->initialize(timeCut, timeMin, timeMax, gradientCut, enableEnergyGradientCut, seedEnergy, minCellEnergy);
88+
mClusterizer->setGeometry(geometry);
89+
LOG(DEBUG) << "Done with clusterizer. Setup cluster factory.";
90+
// Setup cluster factory.
91+
mClusterFactory = decltype(mClusterFactory)(new o2::emcal::ClusterFactory<o2::emcal::Cell>());
92+
LOG(DEBUG) << "Completed init!";
93+
94+
// Setup QA hists.
95+
hCellE.setObject(new TH1F("hCellE", "hCellE", 200, 0.0, 100));
96+
hCellTowerID.setObject(new TH1I("hCellTowerID", "hCellTowerID", 20000, 0, 20000));
97+
hCellEtaPhi.setObject(new TH2F("hCellEtaPhi", "hCellEtaPhi", 160, -0.8, 0.8, 72, 0, 2 * 3.14159));
98+
// NOTE: Reversed column and row because it's more natural for presentatin.
99+
hCellRowCol.setObject(new TH2I("hCellRowCol", "hCellRowCol;Column;Row", 97, 0, 97, 600, 0, 600));
100+
hClusterE.setObject(new TH1F("hClusterE", "hClusterE", 200, 0.0, 100));
101+
hClusterEtaPhi.setObject(new TH2F("hClusterEtaPhi", "hClusterEtaPhi", 160, -0.8, 0.8, 72, 0, 2 * 3.14159));
102+
}
103+
104+
//void process(aod::Collision const& collision, soa::Filtered<aod::Tracks> const& fullTracks, aod::Calos const& cells)
105+
//void process(aod::Collision const& collision, aod::Tracks const& tracks, aod::Calos const& cells)
106+
//void process(aod::BCs const& bcs, aod::Collision const& collision, aod::Calos const& cells)
107+
// Appears to need the BC to be accessed to be available in the collision table...
108+
void process(aod::Collision const& collision, aod::Calos const& cells, aod::BCs const& bcs)
109+
{
110+
LOG(DEBUG) << "Starting process.";
111+
// Convert aod::Calo to o2::emcal::Cell which can be used with the clusterizer.
112+
// In particular, we need to filter only EMCAL cells.
113+
mEmcalCells.clear();
114+
for (auto& cell : cells) {
115+
if (cell.caloType() != selectedCellType || cell.bc() != collision.bc()) {
116+
//LOG(DEBUG) << "Rejected";
117+
continue;
118+
}
119+
//LOG(DEBUG) << "Cell E: " << cell.getEnergy();
120+
//LOG(DEBUG) << "Cell E: " << cell;
121+
122+
mEmcalCells.emplace_back(o2::emcal::Cell(
123+
cell.cellNumber(),
124+
cell.amplitude(),
125+
cell.time(),
126+
o2::emcal::intToChannelType(cell.cellType())));
127+
}
128+
129+
// Cell QA
130+
// For convenience, use the clusterizer stored geometry to get the eta-phi
131+
for (auto& cell : mEmcalCells) {
132+
hCellE->Fill(cell.getEnergy());
133+
hCellTowerID->Fill(cell.getTower());
134+
auto res = mClusterizer->getGeometry()->EtaPhiFromIndex(cell.getTower());
135+
hCellEtaPhi->Fill(std::get<0>(res), std::get<1>(res));
136+
res = mClusterizer->getGeometry()->GlobalRowColFromIndex(cell.getTower());
137+
// NOTE: Reversed column and row because it's more natural for presentatin.
138+
hCellRowCol->Fill(std::get<1>(res), std::get<0>(res));
139+
}
140+
141+
// TODO: Helpful for now, but should be removed.
142+
LOG(DEBUG) << "Converted EMCAL cells";
143+
for (auto& cell : mEmcalCells) {
144+
LOG(DEBUG) << cell.getTower() << ": E: " << cell.getEnergy() << ", time: " << cell.getTimeStamp() << ", type: " << cell.getType();
145+
}
146+
147+
LOG(INFO) << "Converted cells. Contains: " << mEmcalCells.size() << ". Originally " << cells.size() << ". About to run clusterizer.";
148+
149+
// Run the clusterizer
150+
mClusterizer->findClusters(mEmcalCells);
151+
LOG(DEBUG) << "Found clusters.";
152+
auto emcalClusters = mClusterizer->getFoundClusters();
153+
auto emcalClustersInputIndices = mClusterizer->getFoundClustersInputIndices();
154+
LOG(DEBUG) << "Retrieved results. About to setup cluster factory.";
155+
156+
// Convert to analysis clusters.
157+
// First, the cluster factory requires cluster and cell information in order to build the clusters.
158+
mAnalysisClusters.clear();
159+
mClusterFactory->reset();
160+
mClusterFactory->setClustersContainer(*emcalClusters);
161+
mClusterFactory->setCellsContainer(mEmcalCells);
162+
mClusterFactory->setCellsIndicesContainer(*emcalClustersInputIndices);
163+
LOG(DEBUG) << "Cluster factory set up.";
164+
165+
// Convert to analysis clusters.
166+
for (int icl = 0; icl < mClusterFactory->getNumberOfClusters(); icl++) {
167+
auto analysisCluster = mClusterFactory->buildCluster(icl);
168+
mAnalysisClusters.emplace_back(analysisCluster);
169+
}
170+
LOG(DEBUG) << "Converted to analysis clusters.";
171+
172+
// Store the clusters in the table
173+
clusters.reserve(mAnalysisClusters.size());
174+
for (const auto& cluster : mAnalysisClusters) {
175+
// Determine the cluster eta, phi, correcting for the vertex position.
176+
auto pos = cluster.getGlobalPosition();
177+
pos = pos - math_utils::Point3D<float>{collision.posX(), collision.posY(), collision.posZ()};
178+
// Normalize the vector and rescale by energy.
179+
pos /= (cluster.E() / std::sqrt(pos.Mag2()));
180+
181+
// We have our necessary properties. Now we store outputs
182+
//LOG(DEBUG) << "Cluster E: " << cluster.E();
183+
clusters(collision, cluster.E(), pos.Eta(), pos.Phi(), cluster.getM02());
184+
//if (cluster.E() < 0.300) {
185+
// continue;
186+
//}
187+
hClusterE->Fill(cluster.E());
188+
hClusterEtaPhi->Fill(pos.Eta(), pos.Phi());
189+
}
190+
LOG(DEBUG) << "Done with process.";
191+
}
192+
};
193+
194+
WorkflowSpec defineDataProcessing(ConfigContext const&)
195+
{
196+
return WorkflowSpec{
197+
adaptAnalysisTask<EmcalCorrectionTask>("emcal-correction-task")};
198+
}

0 commit comments

Comments
 (0)