Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 33 additions & 21 deletions PWGJE/TableProducer/emcalCorrectionTask.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -191,11 +191,11 @@

// Define the cell energy binning
std::vector<double> cellEnergyBins;
for (int i = 0; i < 51; i++)

Check failure on line 194 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
cellEnergyBins.emplace_back(0.1 * (i - 0) + 0.0); // from 0 to 5 GeV/c, every 0.1 GeV
for (int i = 51; i < 76; i++)

Check failure on line 196 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
cellEnergyBins.emplace_back(0.2 * (i - 51) + 5.2); // from 5.2 to 10.0 GeV, every 0.2 GeV
for (int i = 76; i < 166; i++)

Check failure on line 198 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
cellEnergyBins.emplace_back(1. * (i - 76) + 11.); // from 11.0 to 100. GeV, every 1 GeV

// Setup QA hists.
Expand Down Expand Up @@ -342,7 +342,7 @@

// Store the clusters in the table where a matching collision could
// be identified.
fillClusterTable<CollEventSels::filtered_iterator>(col, vertexPos, iClusterizer, cellIndicesBC, indexMapPair, trackGlobalIndex);
fillClusterTable<CollEventSels::filtered_iterator>(col, vertexPos, iClusterizer, cellIndicesBC, &indexMapPair, &trackGlobalIndex);
} else {
mHistManager.fill(HIST("hBCMatchErrors"), 2);
}
Expand Down Expand Up @@ -470,7 +470,7 @@

// Store the clusters in the table where a matching collision could
// be identified.
fillClusterTable<CollEventSels::filtered_iterator>(col, vertexPos, iClusterizer, cellIndicesBC, indexMapPair, trackGlobalIndex);
fillClusterTable<CollEventSels::filtered_iterator>(col, vertexPos, iClusterizer, cellIndicesBC, &indexMapPair, &trackGlobalIndex);
} else {
mHistManager.fill(HIST("hBCMatchErrors"), 2);
}
Expand Down Expand Up @@ -598,7 +598,7 @@
}
PROCESS_SWITCH(EmcalCorrectionTask, processStandalone, "run stand alone analysis", false);

void cellsToCluster(size_t iClusterizer, const gsl::span<o2::emcal::Cell> cellsBC, std::optional<const gsl::span<o2::emcal::CellLabel>> cellLabels = std::nullopt)
void cellsToCluster(size_t iClusterizer, const gsl::span<o2::emcal::Cell> cellsBC, gsl::span<const o2::emcal::CellLabel> cellLabels = {})
{
mClusterizers.at(iClusterizer)->findClusters(cellsBC);

Expand All @@ -614,10 +614,10 @@
mClusterFactories.reset();
// in preparation for future O2 changes
// mClusterFactories.setClusterizerSettings(mClusterDefinitions.at(iClusterizer).minCellEnergy, mClusterDefinitions.at(iClusterizer).timeMin, mClusterDefinitions.at(iClusterizer).timeMax, mClusterDefinitions.at(iClusterizer).recalcShowerShape5x5);
if (cellLabels) {
mClusterFactories.setContainer(*emcalClusters, cellsBC, *emcalClustersInputIndices, cellLabels);
} else {
if (cellLabels.empty()) {
mClusterFactories.setContainer(*emcalClusters, cellsBC, *emcalClustersInputIndices);
} else {
mClusterFactories.setContainer(*emcalClusters, cellsBC, *emcalClustersInputIndices, cellLabels);
}

LOG(debug) << "Cluster factory set up.";
Expand All @@ -634,24 +634,33 @@
}

template <typename Collision>
void fillClusterTable(Collision const& col, math_utils::Point3D<float> const& vertexPos, size_t iClusterizer, const gsl::span<int64_t> cellIndicesBC, std::optional<std::tuple<std::vector<std::vector<int>>, std::vector<std::vector<int>>>> const& indexMapPair = std::nullopt, std::optional<std::vector<int64_t>> const& trackGlobalIndex = std::nullopt)
void fillClusterTable(Collision const& col, math_utils::Point3D<float> const& vertexPos, size_t iClusterizer, const gsl::span<int64_t> cellIndicesBC, const std::tuple<std::vector<std::vector<int>>, std::vector<std::vector<int>>>* indexMapPair = nullptr, const std::vector<int64_t>* trackGlobalIndex = nullptr)
{
// average number of cells per cluster, only used the reseve a reasonable amount for the clustercells table
const size_t NAvgNcells = 3;

Check failure on line 640 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/function-variable]

Use lowerCamelCase for names of functions and variables.
// we found a collision, put the clusters into the none ambiguous table
clusters.reserve(mAnalysisClusters.size());
if (mClusterLabels.size() > 0) {
if (!mClusterLabels.empty()) {
mcclusters.reserve(mClusterLabels.size());
}
clustercells.reserve(mAnalysisClusters.size() * NAvgNcells);

// get the clusterType once
const auto clusterType = static_cast<int>(mClusterDefinitions[iClusterizer]);

int cellindex = -1;
unsigned int iCluster = 0;
float energy = 0.f;
for (const auto& cluster : mAnalysisClusters) {
energy = cluster.E();
// Determine the cluster eta, phi, correcting for the vertex position.
auto pos = cluster.getGlobalPosition();
pos = pos - vertexPos;
// Normalize the vector and rescale by energy.
pos *= (cluster.E() / std::sqrt(pos.Mag2()));
pos *= (energy / std::sqrt(pos.Mag2()));

// Correct for nonlinear behaviour
float nonlinCorrEnergy = cluster.E();
float nonlinCorrEnergy = energy;
if (!disableNonLin) {
try {
nonlinCorrEnergy = mNonlinearityHandler.getCorrectedClusterEnergy(cluster);
Expand All @@ -662,27 +671,26 @@

// save to table
LOG(debug) << "Writing cluster definition "
<< static_cast<int>(mClusterDefinitions.at(iClusterizer))
<< clusterType
<< " to table.";
mHistManager.fill(HIST("hClusterType"), 1);
clusters(col, cluster.getID(), nonlinCorrEnergy, cluster.getCoreEnergy(), cluster.E(),
clusters(col, cluster.getID(), nonlinCorrEnergy, cluster.getCoreEnergy(), energy,
pos.Eta(), TVector2::Phi_0_2pi(pos.Phi()), cluster.getM02(),
cluster.getM20(), cluster.getNCells(),
cluster.getClusterTime(), cluster.getIsExotic(),
cluster.getDistanceToBadChannel(), cluster.getNExMax(),
static_cast<int>(mClusterDefinitions.at(iClusterizer)));
if (mClusterLabels.size() > 0) {
clusterType);
if (!mClusterLabels.empty()) {
mcclusters(mClusterLabels[iCluster].getLabels(), mClusterLabels[iCluster].getEnergyFractions());
}
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(), cellIndicesBC[cellindex]);
} // end of cells of cluser loop
// fill histograms
mHistManager.fill(HIST("hClusterE"), cluster.E());
mHistManager.fill(HIST("hClusterE"), energy);
mHistManager.fill(HIST("hClusterNLM"), cluster.getNExMax());
mHistManager.fill(HIST("hClusterTime"), cluster.getClusterTime());
mHistManager.fill(HIST("hClusterEtaPhi"), pos.Eta(), TVector2::Phi_0_2pi(pos.Phi()));
Expand All @@ -701,20 +709,25 @@
template <typename BC>
void fillAmbigousClusterTable(BC const& bc, size_t iClusterizer, const gsl::span<int64_t> 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;
int cellindex = -1;
clustersAmbiguous.reserve(mAnalysisClusters.size());
if (mClusterLabels.size() > 0) {
mcclustersAmbiguous.reserve(mClusterLabels.size());
}
clustercellsambiguous.reserve(mAnalysisClusters.size() * NAvgNcells);
unsigned int iCluster = 0;
float energy = 0.f;
for (const auto& cluster : mAnalysisClusters) {
energy = cluster.E();
auto pos = cluster.getGlobalPosition();
pos = pos - math_utils::Point3D<float>{0., 0., 0.};
// Normalize the vector and rescale by energy.
pos *= (cluster.E() / std::sqrt(pos.Mag2()));
pos *= (energy / std::sqrt(pos.Mag2()));

// Correct for nonlinear behaviour
float nonlinCorrEnergy = cluster.E();
float nonlinCorrEnergy = energy;
try {
nonlinCorrEnergy = mNonlinearityHandler.getCorrectedClusterEnergy(cluster);
} catch (o2::emcal::NonlinearityHandler::UninitException& e) {
Expand All @@ -723,22 +736,21 @@

// We have our necessary properties. Now we store outputs

// LOG(debug) << "Cluster E: " << cluster.E();
// LOG(debug) << "Cluster E: " << energy;
if (!hasCollision) {
mHistManager.fill(HIST("hClusterType"), 0);
} else {
mHistManager.fill(HIST("hClusterType"), 2);
}
clustersAmbiguous(
bc, cluster.getID(), nonlinCorrEnergy, cluster.getCoreEnergy(), cluster.E(),
bc, cluster.getID(), nonlinCorrEnergy, cluster.getCoreEnergy(), energy,
pos.Eta(), TVector2::Phi_0_2pi(pos.Phi()), cluster.getM02(),
cluster.getM20(), cluster.getNCells(), cluster.getClusterTime(),
cluster.getIsExotic(), cluster.getDistanceToBadChannel(),
cluster.getNExMax(), static_cast<int>(mClusterDefinitions.at(iClusterizer)));
if (mClusterLabels.size() > 0) {
mcclustersAmbiguous(mClusterLabels[iCluster].getLabels(), mClusterLabels[iCluster].getEnergyFractions());
}
clustercellsambiguous.reserve(cluster.getNCells());
for (int ncell = 0; ncell < cluster.getNCells(); ncell++) {
cellindex = cluster.getCellIndex(ncell);
clustercellsambiguous(clustersAmbiguous.lastIndex(),
Expand Down Expand Up @@ -796,7 +808,7 @@
}
// Tracks that do not point to the EMCal/DCal/PHOS get default values of -999
// This way we can cut out tracks that do not point to the EMCal+DCal
if (track.trackEtaEmcal() < -900 || track.trackPhiEmcal() < -900) {

Check failure on line 811 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
continue;
}
if (trackMinPt > 0 && track.pt() < trackMinPt) {
Expand Down Expand Up @@ -862,7 +874,7 @@
return cellAbsScaleFactors.value[mClusterizers.at(0)->getGeometry()->GetSMType(iSM)];

// Apply cell scale based on columns to accoutn for material of TRD structures
} else if (applyCellAbsScale == 2) {

Check failure on line 877 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
auto res = mClusterizers.at(0)->getGeometry()->GlobalRowColFromIndex(cellID);
return cellAbsScaleFactors.value[std::get<1>(res)];
} else {
Expand All @@ -888,7 +900,7 @@
timeshift = -std::sqrt(215.f + timeCol * timeCol); // 215 is 14.67ns^2 (time it takes to get the cell at eta = 0)

// Also smear the time to account for the broader time resolution in data than in MC
if (cellEnergy < 0.3) // Cells with tless than 300 MeV cannot be the leading cell in the cluster, so their time does not require precise calibration

Check failure on line 903 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
timesmear = 0.; // They will therefore not be smeared and only get their shift
else if (cellType == emcal::ChannelType_t::HIGH_GAIN) // High gain cells -> Low energies
timesmear = normalgaus(rdgen) * (1.6 + 9.5 * std::exp(-3. * cellEnergy)); // Parameters extracted from LHC24f3b & LHC22o (pp), but also usable for other periods
Expand All @@ -896,15 +908,15 @@
timesmear = normalgaus(rdgen) * (5.0); // Parameters extracted from LHC24g4 & LHC24aj (pp), but also usable for other periods

} else { // ---> Data
if (cellEnergy < 0.3) { // Cells with tless than 300 MeV cannot be the leading cell in the cluster, so their time does not require precise calibration

Check failure on line 911 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
timeshift = 0.; // In data they will not be shifted (they are close to 0 anyways)
} else if (cellType == emcal::ChannelType_t::HIGH_GAIN) { // High gain cells -> Low energies
if (cellEnergy < 4.) // Low energy regime

Check failure on line 914 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
timeshift = 0.8 * std::log(2.7 * cellEnergy); // Parameters extracted from LHC22o (pp), but also usable for other periods
else // Medium energy regime
timeshift = 1.5 * std::log(0.9 * cellEnergy); // Parameters extracted from LHC22o (pp), but also usable for other periods
} else if (cellType == emcal::ChannelType_t::LOW_GAIN) { // Low gain cells -> High energies
if (cellEnergy < 30.) // High energy regime

Check failure on line 919 in PWGJE/TableProducer/emcalCorrectionTask.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number]

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
timeshift = 1.9 * std::log(0.09 * cellEnergy); // Parameters extracted from LHC24aj (pp), but also usable for other periods
else // Very high energy regime
timeshift = 1.9; // Parameters extracted from LHC24aj (pp), but also usable for other periods
Expand Down
Loading