diff --git a/PWGJE/Core/CMakeLists.txt b/PWGJE/Core/CMakeLists.txt index 95895ddc442..7eb4bc8ea97 100644 --- a/PWGJE/Core/CMakeLists.txt +++ b/PWGJE/Core/CMakeLists.txt @@ -14,7 +14,8 @@ o2physics_add_library(PWGJECore SOURCES FastJetUtilities.cxx JetFinder.cxx JetBkgSubUtils.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib ONNXRuntime::ONNXRuntime) + emcalCrossTalkEmulation.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib ONNXRuntime::ONNXRuntime O2::EMCALBase O2::EMCALReconstruction) o2physics_target_root_dictionary(PWGJECore HEADERS JetFinder.h @@ -23,5 +24,6 @@ o2physics_target_root_dictionary(PWGJECore JetTaggingUtilities.h JetBkgSubUtils.h JetDerivedDataUtilities.h + emcalCrossTalkEmulation.h LINKDEF PWGJECoreLinkDef.h) endif() diff --git a/PWGJE/Core/emcalCrossTalkEmulation.cxx b/PWGJE/Core/emcalCrossTalkEmulation.cxx new file mode 100644 index 00000000000..9fe9b679461 --- /dev/null +++ b/PWGJE/Core/emcalCrossTalkEmulation.cxx @@ -0,0 +1,604 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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. + +/// \file emcalCrossTalkEmulation.cxx +/// \brief emulation of emcal cross talk for simulations +/// \author Marvin Hemmer , Goethe-University + +#include "emcalCrossTalkEmulation.h" + +#include +#include +#include +#include +#include +#include +#include + +#include // std::find_if +#include +#include // size_t +#include // std::abs +#include // setw +#include // left and right +#include +#include +#include +#include +#include +// #include "Framework/OutputObjHeader.h" + +// #include "Common/CCDB/EventSelectionParams.h" +#include +#include + +#include + +using namespace o2; +using namespace o2::emccrosstalk; +using namespace o2::framework; + +template +auto printArray(std::array const& arr) +{ + std::stringstream ss; + ss << "\n[SM0: " << arr[0]; + for (auto i = 1u; i < N; ++i) { + ss << ", SM" << i << ": " << arr[i]; + } + ss << "]"; + return ss.str(); +} + +template +auto printMatrix(Array2D const& m) +{ + std::stringstream ss; + // Print column headers + ss << std::endl + << std::setw(6) << " " << std::setw(10) << "value1" + << std::setw(10) << "value2" + << std::setw(10) << "value3" + << std::setw(10) << "value4" << std::endl; + + // Print rows with SM labels + for (size_t i = 0; i < m.rows; ++i) { + ss << "SM" << std::left << std::setw(3) << i; // e.g., SM0, SM1... + for (size_t j = 0; j < m.cols; ++j) { + ss << std::right << std::setw(10) << m(i, j); + } + ss << std::endl; + } + + return ss.str(); +} + +void init2DElement(Array2D& matrix, const Array2D& config, const char* name) +{ + int rows = config.rows; + int cols = config.cols; + + if (rows == 0 && cols == 0) { + LOG(info) << name << " has size 0 x 0, so it is disabled!"; + } else if (cols != NNeighbourCases || (rows != 1 && rows != NSM)) { + LOG(error) << name << " must have 4 columns and either 1 or 20 rows!"; + } else { + for (int sm = 0; sm < NSM; ++sm) { + const int row = (rows == 1) ? 0 : sm; + + for (int i = 0; i < cols; ++i) { + matrix[sm][i] = config(row, i); + } + } + } +} + +void init1DElement(std::array& arr, const std::vector& config, const char* name) +{ + size_t confSize = config.size(); + if (confSize == 0) { + LOG(info) << name << " has size 0, so it is disabled!"; + } else if (config.size() != 1 && confSize != NSM) { + LOG(error) << name << " must have either size 1 or 20!"; + } else { + for (int sm = 0; sm < NSM; ++sm) { + const int row = (confSize == 1) ? 0 : sm; + arr[sm] = config[row]; + } + } +} + +o2::framework::AxisSpec axisEnergy = {7000, 0.f, 70.f, "#it{E}_{cell} (GeV)"}; +// For each of the following configurables we will use: +// empty vector == disabled +// vector of size 4 == same for all SM +// vector of vectors with size nSM * 4 == each SM has its own setting +// the 4 values (0-3) correspond to in relative [row,col]: 0: [+-1,0], 1: [+-1,+or-1], 2: [0,+or-1], 3: [+-2, 0 AND +or-1] +// +---+---+-----+---+---+--+--+--+ +// | 3 | 0 | Hit | 0 | 3 | | | | +// +---+---+-----+---+---+--+--+--+ +// | 3 | 1 | 2 | 1 | 3 | | | | +// +---+---+-----+---+---+--+--+--+ + +void EMCCrossTalk::initObjects(const EmcCrossTalkConf& config) +{ + const int run3RunNumber = 223409; + mGeometry = o2::emcal::Geometry::GetInstanceFromRunNumber(run3RunNumber); + if (!mGeometry) { + LOG(error) << "Failure accessing mGeometry"; + } + + // first set the simple run time variables + mTCardCorrClusEnerConserv = config.conserveEnergy.value; + mRandomizeTCard = config.randomizeTCardInducedEnergy.value; + mTCardCorrMinAmp = config.inducedTCardMinimumCellEnergy.value; + mTCardCorrMinInduced = config.inducedTCardMinimum.value; + mTCardCorrMaxInducedELeak = config.inducedTCardMaximumELeak.value; + mTCardCorrMaxInduced = config.inducedTCardMaximum.value; + + // 2nd define the NSM x NNeighbourCases matrices + mTCardCorrInduceEner = Array2D(std::vector(NSM * 4, 0.f), NSM, 4); + mTCardCorrInduceEnerFrac = Array2D(std::vector(NSM * 4, 0.f), NSM, 4); + mTCardCorrInduceEnerFracP1 = Array2D(std::vector(NSM * 4, 0.f), NSM, 4); + mTCardCorrInduceEnerFracWidth = Array2D(std::vector(NSM * 4, 0.f), NSM, 4); + + // now properly init the NSM x NNeighbourCases matrices + // ------------------------------------------------------------------------ + // mTCardCorrInduceEner + init2DElement(mTCardCorrInduceEner, config.inducedEnergyLossConstant.value, "inducedEnergyLossConstant"); + + // mTCardCorrInduceEnerFrac + init2DElement(mTCardCorrInduceEnerFrac, config.inducedEnergyLossFraction.value, "inducedEnergyLossFraction"); + + // mTCardCorrInduceEnerFracP1 + init2DElement(mTCardCorrInduceEnerFracP1, config.inducedEnergyLossFractionP1.value, "inducedEnergyLossFractionP1"); + + // mTCardCorrInduceEnerFracWidth + init2DElement(mTCardCorrInduceEnerFracWidth, config.inducedEnergyLossFractionWidth.value, "inducedEnergyLossFractionWidth"); + // ------------------------------------------------------------------------ + + init1DElement(mTCardCorrInduceEnerFracMax, config.inducedEnergyLossMaximumFraction.value, "inducedEnergyLossMaximumFraction"); + init1DElement(mTCardCorrInduceEnerFracMin, config.inducedEnergyLossMinimumFraction.value, "inducedEnergyLossMinimumFraction"); + init1DElement(mTCardCorrInduceEnerFracMinCentralEta, config.inducedEnergyLossMinimumFractionCentralEta.value, "inducedEnergyLossMinimumFractionCentralEta"); + init1DElement(mTCardCorrInduceEnerProb, config.inducedEnergyLossProbability.value, "inducedEnergyLossProbability"); + + resetArrays(); + + // Print the full matrices and vectors that will be used: + if (config.printConfiguration.value) { + LOGF(info, "inducedEnergyLossConstant: %s", printMatrix((mTCardCorrInduceEner))); + LOGF(info, "inducedEnergyLossFraction: %s", printMatrix((mTCardCorrInduceEnerFrac))); + LOGF(info, "inducedEnergyLossFractionP1: %s", printMatrix((mTCardCorrInduceEnerFracP1))); + LOGF(info, "inducedEnergyLossFractionWidth: %s", printMatrix((mTCardCorrInduceEnerFracWidth))); + LOGF(info, "inducedEnergyLossMaximumFraction: %s", printArray(mTCardCorrInduceEnerFracMax).c_str()); + LOGF(info, "inducedEnergyLossMinimumFraction: %s", printArray(mTCardCorrInduceEnerFracMin).c_str()); + LOGF(info, "inducedEnergyLossMinimumFractionCentralEta: %s", printArray(mTCardCorrInduceEnerFracMinCentralEta).c_str()); + LOGF(info, "inducedEnergyLossProbability: %s", printArray(mTCardCorrInduceEnerProb).c_str()); + } +} + +void EMCCrossTalk::resetArrays() +{ + for (size_t j = 0; j < NCells; j++) { + mTCardCorrCellsEner[j] = 0.; + mTCardCorrCellsNew[j] = false; + } + + mCellsTmp.clear(); +} + +void EMCCrossTalk::setCells(std::vector& cells, std::vector& cellLabels) +{ + mCells = &cells; + mCellsTmp = cells; // a copy since we will need one vector with the changed energies and one with the original ones + mCellLabels = &cellLabels; +} + +void EMCCrossTalk::calculateInducedEnergyInTCardCell(int absId, int absIdRef, int iSM, float ampRef, int cellCase) +{ + // Check that the cell exists + if (absId < 0) { + return; + } + + // Get the fraction + float frac = mTCardCorrInduceEnerFrac[iSM][cellCase] + ampRef * mTCardCorrInduceEnerFracP1[iSM][cellCase]; + + // Use an absolute minimum and maximum fraction if calculated one is out of range + if (frac < mTCardCorrInduceEnerFracMin[iSM]) { + frac = mTCardCorrInduceEnerFracMin[iSM]; + } else if (frac > mTCardCorrInduceEnerFracMax[iSM]) { + frac = mTCardCorrInduceEnerFracMax[iSM]; + } + + // If active, use different absolute minimum fraction for central eta, exclude DCal 2/3 SM + if (mTCardCorrInduceEnerFracMinCentralEta[iSM] > 0 && (iSM < FirstDCal23SM || iSM > LastDCal23SM)) { + // Odd SM + int ietaMin = 32; + int ietaMax = 47; + + // Even SM + if (iSM % 2) { + ietaMin = 0; + ietaMax = 15; + } + + // First get the SM, col-row of this tower + // int imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; + auto [iSM, iMod, iIphi, iIeta] = mGeometry->GetCellIndex(absId); + auto [iphi, ieta] = mGeometry->GetCellPhiEtaIndexInSModule(iSM, iMod, iIphi, iIeta); + + if (ieta >= ietaMin && ieta <= ietaMax) { + if (frac < mTCardCorrInduceEnerFracMinCentralEta[iSM]) + frac = mTCardCorrInduceEnerFracMinCentralEta[iSM]; + } + } // central eta + + LOGF(debug, "\t fraction %2.3f", frac); + + // Randomize the induced fraction, if requested + if (mRandomizeTCard) { + frac = mRandom.Gaus(frac, mTCardCorrInduceEnerFracWidth[iSM][cellCase]); + LOGF(debug, "\t randomized fraction %2.3f", frac); + } + + // If fraction too small or negative, do nothing else + if (frac < Epsilon) { + return; + } + + // Calculate induced energy + float inducedE = mTCardCorrInduceEner[iSM][cellCase] + ampRef * frac; + + // Check if we induce too much energy, in such case use a constant value + if (mTCardCorrMaxInduced < inducedE) + inducedE = mTCardCorrMaxInduced; + + LOGF(debug, "\t induced E %2.3f", inducedE); + + // Try to find the cell that will get energy induced + float amp = 0.f; + auto itCell = std::find_if((*mCells).begin(), (*mCells).end(), [absId](const o2::emcal::Cell& cell) { + return cell.getTower() == absId; + }); + + if (itCell != (*mCells).end()) { + // We found a cell, so let's get the amplitude of that cell + amp = itCell->getAmplitude(); + } else { + amp = 0.f; // this is a new cell, so the base amp is 0.f + } + + // Check that the induced+amp is large enough to avoid extra linearity effects + // typically of the order of the clusterization cell energy cut + // if inducedTCardMaximumELeak was set to a positive value, then induce the energy as long as its smaller than that value + if ((amp + inducedE) > mTCardCorrMinInduced || inducedE < mTCardCorrMaxInducedELeak) { + mTCardCorrCellsEner[absId] += inducedE; + + // If original energy of cell was null, create new one + if (amp <= Epsilon) { + mTCardCorrCellsNew[absId] = true; + } + } else { + return; + } + + LOGF(debug, "Cell %d is with amplitude %2.3f GeV is inducing %1.3f GeV energy to cell %d which already has %2.3f GeV energy with fraction %1.5f", absIdRef, ampRef, inducedE, absId, amp, frac); + + // Subtract the added energy to main cell, if energy conservation is requested + if (mTCardCorrClusEnerConserv) { + mTCardCorrCellsEner[absIdRef] -= inducedE; + } +} + +void EMCCrossTalk::makeCellTCardCorrelation() +{ + int id = -1; + float amp = -1; + + // Loop on all cells with signal + for (const auto& cell : (*mCells)) { + id = cell.getTower(); + amp = cell.getAmplitude(); + + if (amp <= mTCardCorrMinAmp) { + continue; + } + + // First get the SM, col-row of this tower + auto [iSM, iMod, iIphi, iIeta] = mGeometry->GetCellIndex(id); + auto [iphi, ieta] = mGeometry->GetCellPhiEtaIndexInSModule(iSM, iMod, iIphi, iIeta); + + // Determine randomly if we want to create a correlation for this cell, + // depending the SM number of the cell + if (mTCardCorrInduceEnerProb[iSM] < 1) { + if (mRandom.Uniform(0, 1) > mTCardCorrInduceEnerProb[iSM]) { + continue; + } + } + + LOGF(debug, "Reference cell absId %d, iEta %d, iPhi %d, amp %2.3f", id, ieta, iphi, amp); + + // Get the absId of the cells in the cross and same T-Card + int absIDup = -1; + int absIDdo = -1; + int absIDlr = -1; + int absIDuplr = -1; + int absIDdolr = -1; + + int absIDup2 = -1; + int absIDup2lr = -1; + int absIDdo2 = -1; + int absIDdo2lr = -1; + + // Only 2 columns in the T-Card, +1 for even and -1 for odd with respect reference cell + // Sine we only have full T-Cards, we do not need to make any edge case checks + // There is always either a column (eta direction) below or above + int colShift = +1; + if (ieta % 2) { + colShift = -1; + } + + absIDlr = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi, ieta + colShift); + + // Check if up / down cells from reference cell are not out of SM + // First check if there is space one above + if (iphi < emcal::EMCAL_ROWS - 1) { + absIDup = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi + 1, ieta); + absIDuplr = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi + 1, ieta + colShift); + } + + // 2nd check if there is space one below + if (iphi > 0) { + absIDdo = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi - 1, ieta); + absIDdolr = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi - 1, ieta + colShift); + } + + // 3rd check if there is space two above + if (iphi < emcal::EMCAL_ROWS - 2) { + absIDup2 = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi + 2, ieta); + absIDup2lr = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi + 2, ieta + colShift); + } + + // 4th check if there is space two below + if (iphi > 1) { + absIDdo2 = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi - 2, ieta); + absIDdo2lr = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphi - 2, ieta + colShift); + } + + // Check if those cells are in the same T-Card + int tCard = iphi / 8; + if (tCard != (iphi + 1) / 8) { + absIDup = -1; + absIDuplr = -1; + } + if (tCard != (iphi - 1) / 8) { + absIDdo = -1; + absIDdolr = -1; + } + if (tCard != (iphi + 2) / 8) { + absIDup2 = -1; + absIDup2lr = -1; + } + if (tCard != (iphi - 2) / 8) { + absIDdo2 = -1; + absIDdo2lr = -1; + } + + // Calculate induced energy to T-Card cells + // first check if for the given cell case we actually do induce some energy + if (((std::abs(mTCardCorrInduceEner[iSM][0]) > Epsilon) || (std::abs(mTCardCorrInduceEnerFrac[iSM][0]) > Epsilon)) && (std::abs(mTCardCorrInduceEnerFracP1[iSM][0]) > Epsilon) && (std::abs(mTCardCorrInduceEnerFracWidth[iSM][0]) > Epsilon)) { + if (absIDup >= 0) { + LOGF(debug, "cell up %d:", absIDup); + calculateInducedEnergyInTCardCell(absIDup, id, iSM, amp, 0); + } + if (absIDdo >= 0) { + LOGF(debug, "cell down %d:", absIDdo); + calculateInducedEnergyInTCardCell(absIDdo, id, iSM, amp, 0); + } + } + if (((std::abs(mTCardCorrInduceEner[iSM][1]) > Epsilon) || (std::abs(mTCardCorrInduceEnerFrac[iSM][1]) > Epsilon)) && (std::abs(mTCardCorrInduceEnerFracP1[iSM][1]) > Epsilon) && (std::abs(mTCardCorrInduceEnerFracWidth[iSM][1]) > Epsilon)) { + if (absIDuplr >= 0) { + LOGF(debug, "cell up left-right %d:", absIDuplr); + calculateInducedEnergyInTCardCell(absIDuplr, id, iSM, amp, 1); + } + if (absIDdolr >= 0) { + LOGF(debug, "cell down left-right %d:", absIDdolr); + calculateInducedEnergyInTCardCell(absIDdolr, id, iSM, amp, 1); + } + } + if (((std::abs(mTCardCorrInduceEner[iSM][2]) > Epsilon) || (std::abs(mTCardCorrInduceEnerFrac[iSM][2]) > Epsilon)) && (std::abs(mTCardCorrInduceEnerFracP1[iSM][2]) > Epsilon) && (std::abs(mTCardCorrInduceEnerFracWidth[iSM][2]) > Epsilon)) { + if (absIDlr >= 0) { + LOGF(debug, "cell left-right %d:", absIDlr); + calculateInducedEnergyInTCardCell(absIDlr, id, iSM, amp, 2); + } + } + if (((std::abs(mTCardCorrInduceEner[iSM][3]) > Epsilon) || (std::abs(mTCardCorrInduceEnerFrac[iSM][3]) > Epsilon)) && (std::abs(mTCardCorrInduceEnerFracP1[iSM][3]) > Epsilon) && (std::abs(mTCardCorrInduceEnerFracWidth[iSM][3]) > Epsilon)) { + if (absIDup2 >= 0) { + LOGF(debug, "cell up 2nd row %d:", absIDup2); + calculateInducedEnergyInTCardCell(absIDup2, id, iSM, amp, 3); + } + if (absIDdo2 >= 0) { + LOGF(debug, "cell down 2nd row %d:", absIDdo2); + calculateInducedEnergyInTCardCell(absIDdo2, id, iSM, amp, 3); + } + if (absIDup2lr >= 0) { + LOGF(debug, "cell up left-right 2nd row %d:", absIDup2lr); + calculateInducedEnergyInTCardCell(absIDup2lr, id, iSM, amp, 3); + } + if (absIDdo2lr >= 0) { + LOGF(debug, "cell down left-right 2nd row %d:", absIDdo2lr); + calculateInducedEnergyInTCardCell(absIDdo2lr, id, iSM, amp, 3); + } + } + } // cell loop +} + +void EMCCrossTalk::addInducedEnergiesToExistingCells() +{ + // Add the induced energy to the cells and copy them into a new temporal container + // used in AddInducedEnergiesToNewCells() to refill the default cells list fCaloCells + // Create the data member only once. Done here, not sure where to do this properly in the framework. + + for (auto& cell : mCellsTmp) { // o2-linter: disable=const-ref-in-for-loop (we are changing a value here) + float amp = cell.getAmplitude() + mTCardCorrCellsEner[cell.getTower()]; + // Set new amplitude in new temporal container + cell.setAmplitude(amp); + } +} + +void EMCCrossTalk::addInducedEnergiesToNewCells() +{ + // count how many new cells + size_t nCells = (*mCells).size(); + int nCellsNew = 0; + for (size_t j = 0; j < NCells; j++) { + // Newly created? + if (!mTCardCorrCellsNew[j]) { + continue; + } + // Accept only if at least 10 MeV + if (mTCardCorrCellsEner[j] < MinCellEnergy) { + continue; + } + nCellsNew++; + } + + // reserve more space for new cell entries in original cells and celllabels + (*mCells).reserve(nCells + nCellsNew); + (*mCellLabels).reserve(nCells + nCellsNew); + + // change the amplitude of the original cells using + for (size_t iCell = 0; iCell < mCellsTmp.size(); ++iCell) { + (*mCells)[iCell].setAmplitude(mCellsTmp[iCell].getAmplitude()); + } + + // Add the new cells + int absId = -1; + float amp = -1; + float time = 0; + std::vector mclabel; + + for (size_t j = 0; j < NCells; j++) { + // Newly created? + if (!mTCardCorrCellsNew[j]) { + continue; + } + + // Accept only if at least 10 MeV + if (mTCardCorrCellsEner[j] < MinCellEnergy) { + continue; + } + + // Add new cell + absId = j; + amp = mTCardCorrCellsEner[j]; + time = 615.e-9f; + mclabel = {-1}; + + // Assign as MC label the label of the neighboring cell with highest energy + // within the same T-Card. Follow same approach for time. + // Simplest assumption, not fully correct. + // Still assign 0 as fraction of energy. + + // First get the iphi and ieta of this tower + auto [iSM, iMod, iIphi, iIeta] = mGeometry->GetCellIndex(absId); + auto [iphi, ieta] = mGeometry->GetCellPhiEtaIndexInSModule(iSM, iMod, iIphi, iIeta); + + LOGF(debug, "Trying to add cell %d \t ieta = %d\t iphi = %d\t amplitude = %1.3f", absId, ieta, iphi, amp); + + // Loop on the nearest cells around, check the highest energy one, + // and assign its MC label and the time + float ampMax = 0.f; + for (int ietai = ieta - 1; ietai <= ieta + 1; ++ietai) { + for (int iphii = iphi - 1; iphii <= iphi + 1; ++iphii) { + + // Avoid same cell + if (iphii == iphi && ietai == ieta) { + continue; + } + + // Avoid cells out of SM + if (ietai < 0 || ietai >= emcal::EMCAL_COLS || iphii < 0 || iphii >= emcal::EMCAL_ROWS) { + continue; + } + + int absIDi = mGeometry->GetAbsCellIdFromCellIndexes(iSM, iphii, ietai); + // Try to find the cell that will get energy induced + float ampi = 0.f; + size_t indexInCells = 0; + auto itCell = std::find_if((*mCells).begin(), (*mCells).begin() + nCells, [absIDi](const o2::emcal::Cell& cell) { + return cell.getTower() == absIDi; + }); + + if (itCell != (*mCells).begin() + nCells) { + // We found a cell, so let's get the amplitude of that cell + ampi = itCell->getAmplitude(); + if (ampi <= ampMax) { + continue; // early continue if the new amplitude is not the biggest one + } + indexInCells = std::distance((*mCells).begin(), itCell); + LOGF(debug, "Found cell with index %d", indexInCells); + } else { + continue; + } + + // Remove cells with no energy + if (ampi <= MinCellEnergy) { + continue; + } + + // Only same TCard + if (!std::get<0>(mGeometry->areAbsIDsFromSameTCard(absId, absIDi))) { + continue; + } + + std::vector mclabeli = {(*mCellLabels)[indexInCells].GetLeadingMCLabel()}; + float timei = (*mCells)[indexInCells].getTimeStamp(); + + ampMax = ampi; + mclabel = mclabeli; + time = timei; + } // loop phi + } // loop eta + // End Assign MC label + LOGF(debug, "Final ampMax %1.2f\n", ampMax); + LOGF(debug, "--- End : Added cell ID %d, ieta %d, iphi %d, E %1.3f, time %1.3e, mc label %d\n", absId, ieta, iphi, amp, time, mclabel[0]); + + // Add the new cell + (*mCells).emplace_back(absId, amp, time, o2::emcal::intToChannelType(1)); + (*mCellLabels).emplace_back(std::vector{mclabel[0]}, std::vector{0.f}); + } // loop over cells +} + +bool EMCCrossTalk::run() +{ + // START PROCESSING + // Test if cells present + if ((*mCells).size() == 0) { + LOGF(error, "No EMCAL cells found, exiting EMCCrossTalk::run()!"); + return false; + } + + // CELL CROSSTALK EMULATION + // Compute the induced cell energies by T-Card correlation emulation, ONLY MC + makeCellTCardCorrelation(); + + // Add to existing cells the found induced energies in MakeCellTCardCorrelation() if new signal is larger than 10 MeV. + addInducedEnergiesToExistingCells(); + + // Add new cells with found induced energies in MakeCellTCardCorrelation() if new signal is larger than 10 MeV. + addInducedEnergiesToNewCells(); + + resetArrays(); + + return true; +} diff --git a/PWGJE/Core/emcalCrossTalkEmulation.h b/PWGJE/Core/emcalCrossTalkEmulation.h new file mode 100644 index 00000000000..7d83b95d19e --- /dev/null +++ b/PWGJE/Core/emcalCrossTalkEmulation.h @@ -0,0 +1,210 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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. + +/// \file emcalCrossTalkEmulation.h +/// \brief emulation of emcal cross talk for simulations +/// \author Marvin Hemmer , Goethe-University + +#ifndef PWGJE_CORE_EMCALCROSSTALKEMULATION_H_ +#define PWGJE_CORE_EMCALCROSSTALKEMULATION_H_ + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace o2::emccrosstalk +{ +// cell types for enegery induction +enum InductionCellType { + UpDown = 0, + UpDownLeftRight, + LeftRight, + Up2Down2, + NInductionCellType +}; + +// default values for cross talk emulation +// small value to use for comarison equal to 0, std::abs(x) > epsilon +static constexpr float Epsilon = 1e-6f; + +// default for inducedEnergyLossConstant +// static constexpr float DefaultIELC[1][4] = {{0.02f, 0.02f, 0.02f, 0.f}}; // default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml (but is deactivated per default) +static constexpr float DefaultIELC[1][4] = {{0.f, 0.f, 0.f, 0.f}}; // default from pp 13 TeV + +// default for inducedEnergyLossFraction, default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml same as in https://cds.cern.ch/record/2910556/ +static constexpr float DefaultIELF[20][4] = {{1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.20e-02, 1.20e-02, 1.20e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.20e-02, 1.20e-02, 1.20e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.20e-02, 1.20e-02, 1.20e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}, + {1.20e-02, 1.20e-02, 1.20e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}, + {1.15e-02, 1.15e-02, 1.15e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}, + {0.80e-02, 0.80e-02, 0.80e-02, 0.f}}; + +// default for inducedEnergyLossFractionP1, default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml same as in https://cds.cern.ch/record/2910556/ +static constexpr float DefaultIELFP1[1][4] = {{-1.1e-03, -1.1e-03, -1.1e-03, 0.f}}; + +// default for inducedEnergyLossFractionWidth, default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml same as in https://cds.cern.ch/record/2910556/ +static constexpr float DefaultIELFWidth[1][4] = {{5.0e-03, 5.0e-03, 5.0e-03, 0.f}}; + +// default for inducedEnergyLossMinimumFractionCentralEta, IF someone wants to try it. This is purely to document those test numbers from AliEmcalCorrectionConfiguration.yaml! Default is to not use this! Values from default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml +// std::vector DefaultIELMFCE = {6.8e-3, 7.5e-3, 6.8e-3, 9.0e-3, 6.8e-3, 6.8e-3, 6.8e-3, 9.0e-3, 5.2e-3, 5.2e-3, 7.5e-3, 6.8e-3, 6.8e-3, 6.8e-3, 5.2e-3, 5.2e-3, 6.8e-3, 5.2e-3, 5.2e-3, 5.2e-3}; + +struct EmcCrossTalkConf : o2::framework::ConfigurableGroup { + std::string prefix = "emccrosstalk"; + o2::framework::Configurable enableCrossTalk{"enableCrossTalk", false, "Flag to enable cross talk emulation. This should only ever be used for MC!"}; + o2::framework::Configurable createHistograms{"createHistograms", false, "Flag to enable QA histograms."}; + o2::framework::Configurable printConfiguration{"printConfiguration", true, "Flag to print the configuration after initialization."}; + o2::framework::Configurable conserveEnergy{"conserveEnergy", true, "Flag to enable cluster energy conservation."}; + o2::framework::Configurable randomizeTCardInducedEnergy{"randomizeTCardInducedEnergy", true, "Flag to randomize the energy fraction induced by the TCard."}; + o2::framework::Configurable inducedTCardMinimumCellEnergy{"inducedTCardMinimumCellEnergy", 0.01f, "Minimum cell energy in GeV induced by the TCard."}; + o2::framework::Configurable inducedTCardMaximum{"inducedTCardMaximum", 100.f, "Maximum energy in GeV induced by the TCard."}; + o2::framework::Configurable inducedTCardMinimum{"inducedTCardMinimum", 0.1f, "Minimum energy in GeV induced by the TCard + cell energy, IMPORTANT use the same value as the clusterization cell E threshold or not too far from it."}; + o2::framework::Configurable inducedTCardMaximumELeak{"inducedTCardMaximumELeak", 0.f, "Maximum energy in GeV that is going to be leaked independently of what is set with inducedTCardMinimum."}; + // For each of the following Array2D configurables we will use: + // empty vector == disabled + // vector of size 4 == same for all SM + // vector of vectors with size nSM * 4 == each SM has its own setting + // the 4 values (0-3) correspond to in relative [row,col]: 0: [+-1,0], 1: [+-1,+or-1], 2: [0,+or-1], 3: [+-2, 0 AND +or-1] + // +---+---+-----+---+---+--+--+--+ + // | 3 | 0 | Hit | 0 | 3 | | | | + // +---+---+-----+---+---+--+--+--+ + // | 3 | 1 | 2 | 1 | 3 | | | | + // +---+---+-----+---+---+--+--+--+ + // For the std::vector it is similar, empty vector means not used, single value means one value for all SM and 20 values means specifiyng a value for all SM + o2::framework::Configurable> inducedEnergyLossConstant{"inducedEnergyLossConstant", {DefaultIELC[0], 1, 4}, "Constant energy lost by max energy cell in one of T-Card cells. Empty vector == disabled, size 4 vector == enabled. For information on the exact formatting please check the header file."}; + o2::framework::Configurable> inducedEnergyLossFraction{"inducedEnergyLossFraction", {DefaultIELF[0], 20, 4}, "Fraction of energy lost by max energy cell in one of T-Card cells."}; + o2::framework::Configurable> inducedEnergyLossFractionP1{"inducedEnergyLossFractionP1", {DefaultIELFP1[0], 1, 4}, "Slope parameter of fraction of energy lost by max energy cell in one of T-Card cells."}; + o2::framework::Configurable> inducedEnergyLossFractionWidth{"inducedEnergyLossFractionWidth", {DefaultIELFWidth[0], 1, 4}, "Fraction of energy lost by max energy cell in one of T-Card cells, width of random gaussian."}; + // default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml : + // o2::framework::Configurable> inducedEnergyLossMinimumFraction{"inducedEnergyLossMinimumFraction", {3.5e-3f, 5.0e-3f, 4.5e-3f, 6.0e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f, 6.0e-3f, 3.5e-3f, 3.5e-3f, 5.0e-3f, 5.0e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f, 3.5e-3f}, "Minimum induced energy fraction when linear dependency is set."}; + // value from https://cds.cern.ch/record/2910556/: + o2::framework::Configurable> inducedEnergyLossMinimumFraction{"inducedEnergyLossMinimumFraction", {2.35e-3f, 2.5e-3f, 2.35e-3f, 3.0e-3f, 2.35e-3f, 2.35e-3f, 2.35e-3f, 3.0e-3f, 1.75e-3f, 1.75e-3f, 2.5e-3f, 2.35e-3f, 2.35e-3f, 2.35e-3f, 1.75e-3f, 1.75e-3f, 2.35e-3f, 1.75e-3f, 1.75e-3f, 1.75e-3f}, "Minimum induced energy fraction when linear dependency is set."}; + + o2::framework::Configurable> inducedEnergyLossMinimumFractionCentralEta{"inducedEnergyLossMinimumFractionCentralEta", {}, "Minimum induced energy fraction when linear dependency is set. For |eta| < 0.22, if empty no difference in eta. NOT TUNED for TESTING!"}; + + // default from https://github.com/alisw/AliPhysics/blob/master/PWG/EMCAL/config/AliEmcalCorrectionConfiguration.yaml : + // o2::framework::Configurable> inducedEnergyLossMaximumFraction{"inducedEnergyLossMaximumFraction", {0.018f}, "Maximum induced energy fraction when linear dependency is set."}; + // value from https://cds.cern.ch/record/2910556/: + o2::framework::Configurable> inducedEnergyLossMaximumFraction{"inducedEnergyLossMaximumFraction", {0.016f, 0.016f, 0.016f, 0.018f, 0.016f, 0.016f, 0.016f, 0.018f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f, 0.016f}, "Maximum induced energy fraction when linear dependency is set."}; + o2::framework::Configurable> inducedEnergyLossProbability{"inducedEnergyLossProbability", {1.0f}, "Fraction of times max cell energy correlates with cross cells."}; +}; + +static constexpr int NSM = 20; // Number of Supermodules (12 for EMCal + 8 for DCal) +static constexpr int NCells = 17664; // Number of cells in the EMCal +static constexpr int NNeighbourCases = 4; // 0-same row, diff col, 1-up/down cells left/right col 2-left/righ col, and 2nd row cells +static constexpr int FirstDCal23SM = 12; // index of the first 2/3 DCal SM +static constexpr int LastDCal23SM = 17; // index of the last 2/3 DCal SM +static constexpr float MinCellEnergy = 0.01f; // Minimum energy a new cell needs to be added + +// these labels are for later once labeledArrays work on hyperloop. Currently they sadly only allow fixed size not variable size. +// static const std::vector labelsSM{"SM0/all", "SM1", "SM2", "SM3", "SM4", "SM5", "SM6", "SM7", "SM8", "SM9", "SM10", "SM11", "SM12", "SM13", "SM14", "SM15", "SM16", "SM17", "SM18", "SM19"}; +// static const std::vector labelsCells = {"Up&Down", "Up&Down x Left|Right", "Left|Right", "2Up&Down + 2Up&Down xLeft|Right"}; + +class EMCCrossTalk +{ + + public: + ~EMCCrossTalk() + { + LOG(info) << "Destroying EMCCrossTalk"; + } + + /// \brief Basic init function. + /// \param config configurable group containing the config for the cross talk emulation + void initObjects(const EmcCrossTalkConf& config); + + /// \brief Reset arrays containing information for all possible cells. + /// \details mTCardCorrCellsEner and mTCardCorrCellsNew + void resetArrays(); + + /// \brief Sets the pointer the current vector of cells. + /// \param cells pointer to emcal cells of the current event + /// \param cellLabels pointer to emcal cell labels of the current event + void setCells(std::vector& cells, std::vector& cellLabels); + + /// \brief Main function to call later to perform the full cross talk emulation + /// \return flag if everything went well or not + bool run(); + + /// \brief Recover each cell amplitude and absId and induce energy in cells in cross of the same T-Card + void makeCellTCardCorrelation(); + + /// \brief Add to existing cells the found induced energies in makeCellTCardCorrelation() if new signal is larger than 10 MeV. + /// \details Need to destroy/create the default cells list and do a copy from the old to the new via a temporal array fAODCellsTmp. Not too nice or fast, but it works. + void addInducedEnergiesToExistingCells(); + + /// \brief Add new cells with found induced energies in makeCellTCardCorrelation() if new signal is larger than 10 MeV. + void addInducedEnergiesToNewCells(); + + /// \brief Calculate the induced energy in a cell belonging to thesame T-Card as the reference cell. Used in makeCellTCardCorrelation() + /// \param absId Id number of cell in same T-Card as reference cell + /// \param absIdRef Id number of reference cell + /// \param iSM Supermodule number of cell + /// \param ampRef Amplitude of the reference cell + /// \param cellCase Type of cell with respect reference cell 0: up or down, 1: up or down on the diagonal, 2: left or right, 3: 2nd row up/down both left/right + void calculateInducedEnergyInTCardCell(int absId, int absIdRef, int iSM, float ampRef, int cellCase); + + private: + // T-Card correlation emulation, do on MC + bool mTCardCorrClusEnerConserv; // When making correlation, subtract from the reference cell the induced energy on the neighbour cells + std::array mTCardCorrCellsEner; // Array with induced cell energy in T-Card neighbour cells + std::array mTCardCorrCellsNew; // Array with induced cell energy in T-Card neighbour cells, that before had no signal + + o2::framework::Array2D mTCardCorrInduceEner; // Induced energy loss gauss constant on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ col, and 2nd row cells, param 0 + o2::framework::Array2D mTCardCorrInduceEnerFrac; // Induced energy loss gauss fraction param0 on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ col, and 2nd row cells, param 0 + o2::framework::Array2D mTCardCorrInduceEnerFracP1; // Induced energy loss gauss fraction param1 on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ col, and 2nd row cells, param1 + o2::framework::Array2D mTCardCorrInduceEnerFracWidth; // Induced energy loss gauss witdth on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ col, and 2nd row cells + std::array mTCardCorrInduceEnerFracMax; // In case fTCardCorrInduceEnerFracP1 is non null, restrict the maximum fraction of induced energy per SM + std::array mTCardCorrInduceEnerFracMin; // In case fTCardCorrInduceEnerFracP1 is non null, restrict the minimum fraction of induced energy per SM + std::array mTCardCorrInduceEnerFracMinCentralEta; // In case fTCardCorrInduceEnerFracP1 is non null, restrict the minimum fraction of induced energy per SM. Different at central |eta| < 0.22 + std::array mTCardCorrInduceEnerProb; // Probability to induce energy loss per SM + + TRandom3 mRandom; // Random generator + bool mRandomizeTCard; // Use random induced energy + + float mTCardCorrMinAmp; // Minimum cell energy to induce signal on adjacent cells + float mTCardCorrMinInduced; // Minimum induced energy signal on adjacent cells, sum of induced plus original energy, use same as cell energy clusterization cut + float mTCardCorrMaxInducedELeak; // Maximum value of induced energy signal that is always leaked, ~5-10 MeV + float mTCardCorrMaxInduced; // Maximum induced energy signal on adjacent cells + + std::vector* mCells = nullptr; // Pointer to the original cells of the current event + std::vector* mCellLabels = nullptr; // Pointer to the original cell labels of the current event + std::vector mCellsTmp; // Temporal vector of cells (copy) + + o2::emcal::Geometry* mGeometry; // EMCal geometry +}; + +} // namespace o2::emccrosstalk + +#endif // PWGJE_CORE_EMCALCROSSTALKEMULATION_H_ diff --git a/PWGJE/TableProducer/CMakeLists.txt b/PWGJE/TableProducer/CMakeLists.txt index 31e50dab3cf..4644106c740 100644 --- a/PWGJE/TableProducer/CMakeLists.txt +++ b/PWGJE/TableProducer/CMakeLists.txt @@ -93,7 +93,7 @@ endif() o2physics_add_dpl_workflow(emcal-correction-task SOURCES emcalCorrectionTask.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2::EMCALBase O2::EMCALReconstruction O2::EMCALCalibration + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2::EMCALBase O2::EMCALReconstruction O2::EMCALCalibration O2Physics::PWGJECore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(emcal-matchedtracks-writer diff --git a/PWGJE/TableProducer/emcalCorrectionTask.cxx b/PWGJE/TableProducer/emcalCorrectionTask.cxx index b081a9ac5b7..f3c12185e62 100644 --- a/PWGJE/TableProducer/emcalCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalCorrectionTask.cxx @@ -19,6 +19,7 @@ /// #include "PWGJE/Core/JetUtilities.h" +#include "PWGJE/Core/emcalCrossTalkEmulation.h" #include "PWGJE/DataModel/EMCALClusterDefinition.h" #include "PWGJE/DataModel/EMCALClusters.h" #include "PWGJE/DataModel/EMCALMatchedCollisions.h" @@ -26,38 +27,42 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsEMCAL/AnalysisCluster.h" -#include "DataFormatsEMCAL/Cell.h" -#include "DataFormatsEMCAL/CellLabel.h" -#include "DataFormatsEMCAL/Constants.h" -#include "DetectorsBase/GeometryManager.h" -#include "EMCALBase/ClusterFactory.h" -#include "EMCALBase/Geometry.h" -#include "EMCALBase/NonlinearityHandler.h" -#include "EMCALCalib/GainCalibrationFactors.h" -#include "EMCALCalibration/EMCALTempCalibExtractor.h" -#include "EMCALReconstruction/Clusterizer.h" -#include "Framework/ASoA.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" +#include +#include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include +#include #include +#include #include #include #include -#include +#include #include -#include "TVector2.h" #include +#include + +#include #include #include #include #include +#include #include #include #include @@ -70,6 +75,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::emccrosstalk; using MyGlobTracks = o2::soa::Join; using BcEvSels = o2::soa::Join; using CollEventSels = o2::soa::Join; @@ -122,6 +128,12 @@ struct EmcalCorrectionTask { Configurable mcCellEnergyResolutionBroadening{"mcCellEnergyResolutionBroadening", 0., "Relative widening of the MC cell energy resolution. 0 for no widening, 0.1 for 10% widening, etc. Only applied to MC."}; Configurable applyGainCalibShift{"applyGainCalibShift", false, "Apply shift for cell gain calibration to use values before cell format change (Sept. 2023)"}; + // cross talk emulation configs + EmcCrossTalkConf emcCrossTalkConf; + + // cross talk emulation class for handling the cross talk + emccrosstalk::EMCCrossTalk emcCrossTalk; + // Require EMCAL cells (CALO type 1) Filter emccellfilter = aod::calo::caloType == selectedCellType; @@ -162,6 +174,8 @@ struct EmcalCorrectionTask { // Current run number int runNumber{0}; + static constexpr float TrackNotOnEMCal = -900.f; + void init(InitContext const&) { LOG(debug) << "Start init!"; @@ -239,12 +253,15 @@ struct EmcalCorrectionTask { // Define the cell energy binning std::vector cellEnergyBins; - for (int i = 0; i < 51; i++) + for (int i = 0; i < 51; i++) { // o2-linter: disable=magic-number (just numbers for binning) 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++) + } + for (int i = 51; i < 76; i++) { // o2-linter: disable=magic-number (just numbers for binning) 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++) + } + for (int i = 76; i < 166; i++) { // o2-linter: disable=magic-number (just numbers for binning) cellEnergyBins.emplace_back(1. * (i - 76) + 11.); // from 11.0 to 100. GeV, every 1 GeV + } // Setup QA hists. // NOTE: This is not comprehensive. @@ -257,7 +274,8 @@ struct EmcalCorrectionTask { fCrossAxis{100, 0., 1., "F_{+}"}, 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}"}; + nCellAxis{60, -0.5, 59.5, "#it{n}_{cells}"}, + energyDenseAxis = {7000, 0.f, 70.f, "#it{E}_{cell} (GeV)"}; mHistManager.add("hCellE", "hCellE", O2HistType::kTH1D, {energyAxis}); mHistManager.add("hCellTowerID", "hCellTowerID", O2HistType::kTH1D, {{20000, 0, 20000}}); mHistManager.add("hCellEtaPhi", "hCellEtaPhi", O2HistType::kTH2F, {etaAxis, phiAxis}); @@ -310,6 +328,14 @@ struct EmcalCorrectionTask { mHistManager.add("hClusterFCrossSigmaShortE", "hClusterFCrossSigmaShortE", O2HistType::kTH3F, {energyAxis, fCrossAxis, sigmaShortAxis}); } + if (isMC.value && emcCrossTalkConf.enableCrossTalk.value) { + emcCrossTalk.initObjects(emcCrossTalkConf); + if (emcCrossTalkConf.createHistograms.value) { + mHistManager.add("hCellEnergyDistBefore", "Cell energy before cross-talk emulation", {o2::framework::HistType::kTH1D, {energyDenseAxis}}); + mHistManager.add("hCellEnergyDistAfter", "Cell energy after cross-talk emulation", {o2::framework::HistType::kTH1D, {energyDenseAxis}}); + } + } + // For some runs, LG cells require an extra time shift of 2 * 8.8ns due to problems in the time calibration // Affected run ranges (inclusive) are initialised here (min,max) mExtraTimeShiftRunRanges.emplace_back(535365, 535645); // LHC23g-LHC23h @@ -500,6 +526,13 @@ struct EmcalCorrectionTask { mHistManager.fill(HIST("hContributors"), cell.mcParticle_as().size()); auto cellParticles = cell.mcParticle_as(); for (const auto& cellparticle : cellParticles) { + const auto& ids = cell.mcParticleIds(); + const auto& amps = cell.amplitudeA(); + + if (ids.empty() || amps.empty()) { + LOGF(warning, "Skipping cell with empty MC info: absId=%d", cell.cellNumber()); + continue; + } mHistManager.fill(HIST("hMCParticleEnergy"), cellparticle.e()); } auto amplitude = cell.amplitude(); @@ -517,7 +550,39 @@ struct EmcalCorrectionTask { cell.time() + getCellTimeShift(cell.cellNumber(), amplitude, o2::emcal::intToChannelType(cell.cellType()), runNumber), o2::emcal::intToChannelType(cell.cellType())); cellIndicesBC.emplace_back(cell.globalIndex()); - cellLabels.emplace_back(cell.mcParticleIds(), cell.amplitudeA()); + cellLabels.emplace_back(std::vector{cell.mcParticleIds().begin(), cell.mcParticleIds().end()}, std::vector{cell.amplitudeA().begin(), cell.amplitudeA().end()}); + } + if (isMC.value && emcCrossTalkConf.enableCrossTalk.value) { + if (emcCrossTalkConf.createHistograms.value) { + for (const auto& cell : cellsBC) { + mHistManager.fill(HIST("hCellEnergyDistBefore"), cell.getAmplitude()); + } + } + emcCrossTalk.setCells(cellsBC, cellLabels); + bool isOkCrossTalk = emcCrossTalk.run(); + if (!isOkCrossTalk) { + LOG(info) << "Cross talk emulation failed!"; + } else { + // When we get new cells we also need to add additional entries into cellIndicesBC. + // Adding -1 and later when filling the clusterID<->cellID table skip all cases where this is -1 + if (cellIndicesBC.size() < cellsBC.size()) { + cellIndicesBC.reserve(cellsBC.size()); + for (size_t iMissing = 0; iMissing < (cellsBC.size() - cellIndicesBC.size()); ++iMissing) { + cellIndicesBC.emplace_back(-1); + } + } + if (emcCrossTalkConf.createHistograms.value) { + for (const auto& cell : cellsBC) { + mHistManager.fill(HIST("hCellEnergyDistAfter"), cell.getAmplitude()); + } + } + } // cross talk emulation was okay + } // if (isMC.value && emcCrossTalkConf.enableCrossTalk.value) + // shaper correction has to come AFTER cross talk + for (auto& cell : cellsBC) { // o2-linter: disable=const-ref-in-for-loop (we are changing a value here) + if (cell.getLowGain()) { + cell.setAmplitude(o2::emcal::NonlinearityHandler::evaluateShaperCorrectionCellEnergy(cell.getAmplitude())); + } } LOG(detail) << "Number of cells for BC (CF): " << cellsBC.size(); nCellsProcessed += cellsBC.size(); @@ -786,7 +851,9 @@ struct EmcalCorrectionTask { 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]); + if (cellIndicesBC[cellindex] >= 0) { + clustercells(clusters.lastIndex(), cellIndicesBC[cellindex]); + } } // end of cells of cluser loop // fill histograms mHistManager.fill(HIST("hClusterE"), energy); @@ -907,13 +974,12 @@ struct EmcalCorrectionTask { { int nTrack = 0; for (const auto& track : tracks) { - // TODO only consider tracks in current emcal/dcal acceptanc if (!track.isGlobalTrack()) { // only global tracks continue; } // 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) { + if (track.trackEtaEmcal() < TrackNotOnEMCal || track.trackPhiEmcal() < TrackNotOnEMCal) { continue; } if (trackMinPt > 0 && track.pt() < trackMinPt) {