From 2bf38b0c7fe5988b8fdfd866de3c9d39c89a9509 Mon Sep 17 00:00:00 2001 From: peressounko Date: Fri, 17 Jan 2025 13:18:47 +0300 Subject: [PATCH 1/2] ClusterProducer: Final non-linearity added; SW trigger selection added --- Common/TableProducer/caloClusterProducer.cxx | 352 +++++++++--------- PWGEM/Tasks/CMakeLists.txt | 4 +- PWGEM/Tasks/phosNonlin.cxx | 150 ++++---- PWGEM/Tasks/phosPi0.cxx | 353 ++++++++++++------- 4 files changed, 497 insertions(+), 362 deletions(-) diff --git a/Common/TableProducer/caloClusterProducer.cxx b/Common/TableProducer/caloClusterProducer.cxx index b569cf7a58f..a6afbd901db 100644 --- a/Common/TableProducer/caloClusterProducer.cxx +++ b/Common/TableProducer/caloClusterProducer.cxx @@ -9,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file caloClusterProducer.cxx +/// \brief Produces PHOS clusters from PHOS cells +/// +/// \author Dmitri Peresunko + #include "Framework/ConfigParamSpec.h" #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" @@ -36,7 +41,7 @@ using namespace o2::framework; using namespace o2; -using mcCells = o2::soa::Join; +using McCells = o2::soa::Join; struct caloClusterProducerTask { Produces clucursor; @@ -46,13 +51,13 @@ struct caloClusterProducerTask { Produces cluambmccursor; Configurable isMC{"isMC", 0, "0 - data, 1 - MC"}; - Configurable useCoreE{"coreE", 0, "0 - full energy, 1 - core energy"}; + Configurable useCoreE{"useCoreE", 0, "0 - full energy, 1 - core energy"}; Configurable skipL1phase{"skipL1phase", false, "skip or apply L1phase time correction"}; - Configurable mNonlinType{"nonlinType", 1, "0:no corr, 1: default"}; - Configurable> cpvMinE{"cpvCluMinAmp", {20., 50., 50.}, "minimal CPV cluster amplitude per module"}; - Configurable mBadMapPath{"badmapPath", "PHS/Calib/BadMap", "path to BadMap snapshot"}; - Configurable mCalibPath{"calibPath", "PHS/Calib/CalibParams", "path to Calibration snapshot"}; - Configurable mL1PhasePath{"L1phasePath", "PHS/Calib/L1phase", "path to L1phase snapshot"}; + Configurable nonlinType{"nonlinType", 1, "0:no corr, 1: data, 2: MC"}; + Configurable> cpvMinE{"cpvMinE", {20., 50., 50.}, "minimal CPV cluster amplitude per module"}; + Configurable badMapPath{"badMapPath", "PHS/Calib/BadMap", "path to BadMap snapshot"}; + Configurable calibPath{"calibPath", "PHS/Calib/CalibParams", "path to Calibration snapshot"}; + Configurable l1phasePath{"l1phasePath", "PHS/Calib/L1phase", "path to L1phase snapshot"}; Service ccdb; @@ -71,15 +76,15 @@ struct caloClusterProducerTask { static constexpr int16_t kCpvX = 7; // grid 6 steps along z and 7 along phi as largest match ellips 20x20 cm static constexpr int16_t kCpvZ = 6; static constexpr int16_t kCpvCells = 4 * kCpvX * kCpvZ; // 4 modules - static constexpr float cpvMaxX = 73; // max CPV coordinate phi - static constexpr float cpvMaxZ = 63; // max CPV coordinate z + static constexpr float kCpvMaxX = 73; // max CPV coordinate phi + static constexpr float kCpvMaxZ = 63; // max CPV coordinate z - class trackMatch + class TrackMatch { public: - trackMatch() = default; - trackMatch(float x, float z, int i) : pX(x), pZ(z), indx(i) {} - ~trackMatch() = default; + TrackMatch() = default; + TrackMatch(float x, float z, int i) : pX(x), pZ(z), indx(i) {} + ~TrackMatch() = default; public: float pX = 9999.; // X (phi) track coordinate in PHOS plane @@ -87,11 +92,11 @@ struct caloClusterProducerTask { int indx = -1; // track global index }; - class trackTrigRec + class TrackTrigRec { public: - trackTrigRec() = default; - ~trackTrigRec() = default; + TrackTrigRec() = default; + ~TrackTrigRec() = default; public: int64_t mTR; // BC ref @@ -123,7 +128,7 @@ struct caloClusterProducerTask { } std::map bcMap; int bcId = 0; - for (auto bc : bcs) { + for (auto & bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -131,7 +136,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (auto cl : colls) { + for (const auto & cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -149,11 +154,11 @@ struct caloClusterProducerTask { // Fill output table // calibration may be updated by CCDB fetcher - const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(mBadMapPath, timestamp); - const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(mCalibPath, timestamp); + const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(badMapPath, timestamp); + const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(calibPath, timestamp); if (!isMC && !skipL1phase) { - const std::vector* vec = ccdb->getForTimeStamp>(mL1PhasePath, timestamp); + const std::vector* vec = ccdb->getForTimeStamp>(l1phasePath, timestamp); if (vec) { clusterizerPHOS->setL1phase((*vec)[0]); } else { @@ -182,7 +187,7 @@ struct caloClusterProducerTask { o2::InteractionRecord ir; const int kPHOS = 0; - for (auto& c : cells) { + for (const auto& c : cells) { if (c.caloType() != kPHOS) // PHOS continue; // Fix for bug in trigger digits @@ -216,7 +221,7 @@ struct caloClusterProducerTask { // Find CPV clusters corresponding to PHOS trigger records std::vector> cpvMatchPoints[kCpvCells]; // Number of entries in each cell per TrigRecord - std::vector cpvNMatchPoints; + std::vector cpvNMatchPoints; cpvNMatchPoints.reserve(outputPHOSClusterTrigRecs.size()); int64_t curBC = -1; if (cpvs.begin() != cpvs.end()) { @@ -245,7 +250,7 @@ struct caloClusterProducerTask { if (cpvclu.amplitude() < static_cast>(cpvMinE)[static_cast(cpvclu.moduleNumber()) - 2]) { continue; } - int index = CpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); + int index = cpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); cpvMatchPoints[index].emplace_back(cpvclu.posX(), cpvclu.posZ()); } if (cpvNMatchPoints.size() > 0) { @@ -255,7 +260,7 @@ struct caloClusterProducerTask { } // Fill output - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { int firstClusterInEvent = cluTR.getFirstEntry(); int lastClusterInEvent = firstClusterInEvent + cluTR.getNumberOfObjects(); @@ -293,7 +298,7 @@ struct caloClusterProducerTask { // Correction for the depth of the shower starting point (TDR p 127) const float para = 0.925; const float parb = 6.52; - float depth = para * TMath::Log(e) + parb; + float depth = para * std::log(e) + parb; posX -= posX * depth / 460.; posZ -= (posZ - vtx.Z()) * depth / 460.; @@ -306,51 +311,51 @@ struct caloClusterProducerTask { continue; } - e = Nonlinearity(e); + e = nonlinearity(e); mom.SetMag(e); float cpvdist = 99.; - const float cellSizeX = 2 * cpvMaxX / kCpvX; - const float cellSizeZ = 2 * cpvMaxZ / kCpvZ; + const float cellSizeX = 2 * kCpvMaxX / kCpvX; + const float cellSizeZ = 2 * kCpvMaxZ / kCpvZ; // look 9 CPV regions around PHOS cluster if (mod >= 2 && cpvExist) { // CPV exist in mods 2,3,4 - int phosIndex = CpvMatchIndex(mod, posX, posZ); + int phosIndex = cpvMatchIndex(mod, posX, posZ); std::vector regions; regions.push_back(phosIndex); - if (posX > -cpvMaxX + cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom left + if (posX > -kCpvMaxX + cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom left regions.push_back(phosIndex - kCpvZ - 1); } regions.push_back(phosIndex - kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top left + if (posZ < kCpvMaxZ - cellSizeZ) { // top left regions.push_back(phosIndex - kCpvZ + 1); } } - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom regions.push_back(phosIndex - 1); } - if (posZ < cpvMaxZ - cellSizeZ) { // top + if (posZ < kCpvMaxZ - cellSizeZ) { // top regions.push_back(phosIndex + 1); } - if (posX < cpvMaxX - cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom right + if (posX < kCpvMaxX - cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom right regions.push_back(phosIndex + kCpvZ - 1); } regions.push_back(phosIndex + kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top right + if (posZ < kCpvMaxZ - cellSizeZ) { // top right regions.push_back(phosIndex + kCpvZ + 1); } } - float sigmaX = 1. / TMath::Min(5.2, 1.111 + 0.56 * TMath::Exp(-0.031 * e * e) + 4.8 / TMath::Power(e + 0.61, 3)); // inverse sigma X - float sigmaZ = 1. / TMath::Min(3.3, 1.12 + 0.35 * TMath::Exp(-0.032 * e * e) + 0.75 / TMath::Power(e + 0.24, 3)); // inverse sigma Z + float sigmaX = 1. / std::min(5.2, 1.111 + 0.56 * std::exp(-0.031 * e * e) + 4.8 / std::pow(e + 0.61, 3)); // inverse sigma X + float sigmaZ = 1. / std::min(3.3, 1.12 + 0.35 * std::exp(-0.032 * e * e) + 0.75 / std::pow(e + 0.24, 3)); // inverse sigma Z - for (int indx : regions) { + for (const int & indx : regions) { if (indx >= 0 && indx < kCpvCells) { for (int ii = cpvPoints->mStart[indx]; ii < cpvPoints->mEnd[indx]; ii++) { auto p = cpvMatchPoints[indx][ii]; - float d = pow((p.first - posX) * sigmaX, 2) + pow((p.second - posZ) * sigmaZ, 2); + float d = std::pow((p.first - posX) * sigmaX, 2) + std::pow((p.second - posZ) * sigmaZ, 2); if (d < cpvdist) { cpvdist = d; } @@ -359,7 +364,7 @@ struct caloClusterProducerTask { } } if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + cpvdist = std::sqrt(cpvdist); // was squared } int cpvindex = -2; // -2 no CPV in event if (cpvExist) { @@ -404,7 +409,7 @@ struct caloClusterProducerTask { void processStandaloneMC(o2::aod::BCsWithTimestamps const& bcs, o2::aod::Collisions const& colls, - mcCells& cells, + McCells const& cells, o2::aod::CaloTriggers const&, o2::aod::CPVClusters const& cpvs) { @@ -417,7 +422,7 @@ struct caloClusterProducerTask { } std::map bcMap; int bcId = 0; - for (auto bc : bcs) { + for (auto const & bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -425,7 +430,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (auto cl : colls) { + for (auto const &cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -443,11 +448,11 @@ struct caloClusterProducerTask { // Fill output table // calibration may be updated by CCDB fetcher - const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(mBadMapPath, timestamp); - const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(mCalibPath, timestamp); + const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(badMapPath, timestamp); + const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(calibPath, timestamp); if (!isMC && !skipL1phase) { - const std::vector* vec = ccdb->getForTimeStamp>(mL1PhasePath, timestamp); + const std::vector* vec = ccdb->getForTimeStamp>(l1phasePath, timestamp); if (vec) { clusterizerPHOS->setL1phase((*vec)[0]); } else { @@ -477,7 +482,7 @@ struct caloClusterProducerTask { o2::InteractionRecord ir; const int kPHOS = 0; o2::dataformats::MCTruthContainer cellTruth; - for (auto& c : cells) { + for (const auto& c : cells) { if (c.caloType() != kPHOS) // PHOS continue; // Fix for bug in trigger digits @@ -532,7 +537,7 @@ struct caloClusterProducerTask { // Find CPV clusters corresponding to PHOS trigger records std::vector> cpvMatchPoints[kCpvCells]; // Number of entries in each cell per TrigRecord - std::vector cpvNMatchPoints; + std::vector cpvNMatchPoints; cpvNMatchPoints.reserve(outputPHOSClusterTrigRecs.size()); int64_t curBC = -1; if (cpvs.begin() != cpvs.end()) { @@ -561,7 +566,7 @@ struct caloClusterProducerTask { if (cpvclu.amplitude() < static_cast>(cpvMinE)[static_cast(cpvclu.moduleNumber()) - 2]) { continue; } - int index = CpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); + int index = cpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); cpvMatchPoints[index].emplace_back(cpvclu.posX(), cpvclu.posZ()); } if (cpvNMatchPoints.size() > 0) { @@ -571,7 +576,7 @@ struct caloClusterProducerTask { } // Fill output - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { int firstClusterInEvent = cluTR.getFirstEntry(); int lastClusterInEvent = firstClusterInEvent + cluTR.getNumberOfObjects(); @@ -609,7 +614,7 @@ struct caloClusterProducerTask { // Correction for the depth of the shower starting point (TDR p 127) const float para = 0.925; const float parb = 6.52; - float depth = para * TMath::Log(e) + parb; + float depth = para * std::log(e) + parb; posX -= posX * depth / 460.; posZ -= (posZ - vtx.Z()) * depth / 460.; @@ -622,51 +627,51 @@ struct caloClusterProducerTask { continue; } - e = Nonlinearity(e); + e = nonlinearity(e); mom.SetMag(e); float cpvdist = 99.; - const float cellSizeX = 2 * cpvMaxX / kCpvX; - const float cellSizeZ = 2 * cpvMaxZ / kCpvZ; + const float cellSizeX = 2 * kCpvMaxX / kCpvX; + const float cellSizeZ = 2 * kCpvMaxZ / kCpvZ; // look 9 CPV regions around PHOS cluster if (mod >= 2 && cpvExist) { // CPV exist in mods 2,3,4 - int phosIndex = CpvMatchIndex(mod, posX, posZ); + int phosIndex = cpvMatchIndex(mod, posX, posZ); std::vector regions; regions.push_back(phosIndex); - if (posX > -cpvMaxX + cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom left + if (posX > -kCpvMaxX + cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom left regions.push_back(phosIndex - kCpvZ - 1); } regions.push_back(phosIndex - kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top left + if (posZ < kCpvMaxZ - cellSizeZ) { // top left regions.push_back(phosIndex - kCpvZ + 1); } } - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom regions.push_back(phosIndex - 1); } - if (posZ < cpvMaxZ - cellSizeZ) { // top + if (posZ < kCpvMaxZ - cellSizeZ) { // top regions.push_back(phosIndex + 1); } - if (posX < cpvMaxX - cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom right + if (posX < kCpvMaxX - cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom right regions.push_back(phosIndex + kCpvZ - 1); } regions.push_back(phosIndex + kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top right + if (posZ < kCpvMaxZ - cellSizeZ) { // top right regions.push_back(phosIndex + kCpvZ + 1); } } - float sigmaX = 1. / TMath::Min(5.2, 1.111 + 0.56 * TMath::Exp(-0.031 * e * e) + 4.8 / TMath::Power(e + 0.61, 3)); // inverse sigma X - float sigmaZ = 1. / TMath::Min(3.3, 1.12 + 0.35 * TMath::Exp(-0.032 * e * e) + 0.75 / TMath::Power(e + 0.24, 3)); // inverse sigma Z + float sigmaX = 1. / std::min(5.2, 1.111 + 0.56 * std::exp(-0.031 * e * e) + 4.8 / std::pow(e + 0.61, 3)); // inverse sigma X + float sigmaZ = 1. / std::min(3.3, 1.12 + 0.35 * std::exp(-0.032 * e * e) + 0.75 / std::pow(e + 0.24, 3)); // inverse sigma Z - for (int indx : regions) { + for (const int indx : regions) { if (indx >= 0 && indx < kCpvCells) { for (int ii = cpvPoints->mStart[indx]; ii < cpvPoints->mEnd[indx]; ii++) { auto p = cpvMatchPoints[indx][ii]; - float d = pow((p.first - posX) * sigmaX, 2) + pow((p.second - posZ) * sigmaZ, 2); + float d = std::pow((p.first - posX) * sigmaX, 2) + std::pow((p.second - posZ) * sigmaZ, 2); if (d < cpvdist) { cpvdist = d; } @@ -675,7 +680,7 @@ struct caloClusterProducerTask { } } if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + cpvdist = std::sqrt(cpvdist); // was squared } int cpvindex = -2; // -2 no CPV in event if (cpvExist) { @@ -689,7 +694,7 @@ struct caloClusterProducerTask { mclabels.clear(); mcamplitudes.clear(); gsl::span spDigList = outputTruthCont.getLabels(i); - for (auto cellLab : spDigList) { + for (const auto &cellLab : spDigList) { mclabels.push_back(cellLab.getTrackID()); // Track ID in current event? mcamplitudes.push_back(cellLab.getEdep()); } @@ -760,7 +765,7 @@ struct caloClusterProducerTask { std::map bcMap; int bcId = 0; - for (auto bc : bcs) { + for (const auto &bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -768,7 +773,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (auto cl : colls) { + for (const auto &cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -785,11 +790,11 @@ struct caloClusterProducerTask { // Fill output table // calibration may be updated by CCDB fetcher - const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(mBadMapPath, timestamp); - const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(mCalibPath, timestamp); + const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(badMapPath, timestamp); + const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(calibPath, timestamp); if (!isMC && !skipL1phase) { - const std::vector* vec = ccdb->getForTimeStamp>(mL1PhasePath, timestamp); + const std::vector* vec = ccdb->getForTimeStamp>(l1phasePath, timestamp); if (vec) { clusterizerPHOS->setL1phase((*vec)[0]); } else { @@ -818,7 +823,7 @@ struct caloClusterProducerTask { o2::InteractionRecord ir; const int kPHOS = 0; - for (auto& c : cells) { + for (const auto& c : cells) { if (c.caloType() != kPHOS) // PHOS continue; // Fix for bug in trigger digits @@ -852,7 +857,7 @@ struct caloClusterProducerTask { // Find CPV clusters corresponding to PHOS trigger records std::vector> cpvMatchPoints[kCpvCells]; // Number of entries in each cell per TrigRecord - std::vector cpvNMatchPoints; + std::vector cpvNMatchPoints; cpvNMatchPoints.reserve(outputPHOSClusterTrigRecs.size()); int64_t curBC = -1; @@ -881,7 +886,7 @@ struct caloClusterProducerTask { if (cpvclu.amplitude() < static_cast>(cpvMinE)[static_cast(cpvclu.moduleNumber()) - 2]) { continue; } - int index = CpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); + int index = cpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); cpvMatchPoints[index].emplace_back(cpvclu.posX(), cpvclu.posZ()); } if (cpvNMatchPoints.size()) { @@ -890,9 +895,9 @@ struct caloClusterProducerTask { } } // same for tracks - std::vector trackMatchPoints[kCpvCells]; // tracks hit in grid/cell in PHOS + std::vector trackMatchPoints[kCpvCells]; // tracks hit in grid/cell in PHOS // Number of entries in each cell per TrigRecord - std::vector trackNMatchPoints; + std::vector trackNMatchPoints; trackNMatchPoints.reserve(outputPHOSClusterTrigRecs.size()); curBC = 0; @@ -903,7 +908,7 @@ struct caloClusterProducerTask { } } bool keepBC = false; - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { if (cluTR.getBCData().toLong() == curBC) { keepBC = true; break; @@ -930,7 +935,7 @@ struct caloClusterProducerTask { curBC = track.collision().bc_as().globalBC(); } keepBC = false; - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { if (cluTR.getBCData().toLong() == curBC) { keepBC = true; break; @@ -954,7 +959,7 @@ struct caloClusterProducerTask { float trackX, trackZ; auto trackPar = getTrackPar(track); if (impactOnPHOS(trackPar, track.trackEtaEmcal(), track.trackPhiEmcal(), track.collision().posZ(), module, trackX, trackZ)) { - int index = CpvMatchIndex(module, trackX, trackZ); + int index = cpvMatchIndex(module, trackX, trackZ); trackMatchPoints[index].emplace_back(trackX, trackZ, track.globalIndex()); } } @@ -965,7 +970,7 @@ struct caloClusterProducerTask { } // Fill output tables - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { int firstClusterInEvent = cluTR.getFirstEntry(); int lastClusterInEvent = firstClusterInEvent + cluTR.getNumberOfObjects(); @@ -1012,7 +1017,7 @@ struct caloClusterProducerTask { // Correction for the depth of the shower starting point (TDR p 127) const float para = 0.925; const float parb = 6.52; - float depth = para * TMath::Log(e) + parb; + float depth = para * std::log(e) + parb; posX -= posX * depth / 460.; posZ -= (posZ - vtx.Z()) * depth / 460.; @@ -1025,53 +1030,53 @@ struct caloClusterProducerTask { continue; } - e = Nonlinearity(e); + e = nonlinearity(e); mom.SetMag(e); // CPV and track match - const float cellSizeX = 2 * cpvMaxX / kCpvX; - const float cellSizeZ = 2 * cpvMaxZ / kCpvZ; + const float cellSizeX = 2 * kCpvMaxX / kCpvX; + const float cellSizeZ = 2 * kCpvMaxZ / kCpvZ; // look 9 CPV regions around PHOS cluster - int phosIndex = CpvMatchIndex(mod, posX, posZ); + int phosIndex = cpvMatchIndex(mod, posX, posZ); std::vector regions; regions.push_back(phosIndex); - if (posX > -cpvMaxX + cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom left + if (posX > -kCpvMaxX + cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom left regions.push_back(phosIndex - kCpvZ - 1); } regions.push_back(phosIndex - kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top left + if (posZ < kCpvMaxZ - cellSizeZ) { // top left regions.push_back(phosIndex - kCpvZ + 1); } } - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom regions.push_back(phosIndex - 1); } - if (posZ < cpvMaxZ - cellSizeZ) { // top + if (posZ < kCpvMaxZ - cellSizeZ) { // top regions.push_back(phosIndex + 1); } - if (posX < cpvMaxX - cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom right + if (posX < kCpvMaxX - cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom right regions.push_back(phosIndex + kCpvZ - 1); } regions.push_back(phosIndex + kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top right + if (posZ < kCpvMaxZ - cellSizeZ) { // top right regions.push_back(phosIndex + kCpvZ + 1); } } - float sigmaX = 1. / TMath::Min(5.2, 1.111 + 0.56 * TMath::Exp(-0.031 * e * e) + 4.8 / TMath::Power(e + 0.61, 3)); // inverse sigma X - float sigmaZ = 1. / TMath::Min(3.3, 1.12 + 0.35 * TMath::Exp(-0.032 * e * e) + 0.75 / TMath::Power(e + 0.24, 3)); // inverse sigma Z + float sigmaX = 1. / std::min(5.2, 1.111 + 0.56 * std::exp(-0.031 * e * e) + 4.8 / std::pow(e + 0.61, 3)); // inverse sigma X + float sigmaZ = 1. / std::min(3.3, 1.12 + 0.35 * std::exp(-0.032 * e * e) + 0.75 / std::pow(e + 0.24, 3)); // inverse sigma Z float cpvdist = 99., trackdist = 99.; // float cpvDx = 0., cpvDz = 0.; float trackDx = 9999., trackDz = 9999.; int trackindex = -1; - for (int indx : regions) { + for (const int indx : regions) { if (cpvPoints != cpvNMatchPoints.end()) { if (indx >= 0 && indx < kCpvCells) { for (int ii = cpvPoints->mStart[indx]; ii < cpvPoints->mEnd[indx]; ii++) { auto p = cpvMatchPoints[indx][ii]; - float d = pow((p.first - posX) * sigmaX, 2) + pow((p.second - posZ) * sigmaZ, 2); + float d = std::pow((p.first - posX) * sigmaX, 2) + std::pow((p.second - posZ) * sigmaZ, 2); if (d < cpvdist) { cpvdist = d; } @@ -1083,7 +1088,7 @@ struct caloClusterProducerTask { if (trackPoints != trackNMatchPoints.end()) { for (int ii = trackPoints->mStart[indx]; ii < trackPoints->mEnd[indx]; ii++) { auto pp = trackMatchPoints[indx][ii]; - float d = pow((pp.pX - posX) * sigmaX, 2) + pow((pp.pZ - posZ) * sigmaZ, 2); // TODO different sigma for tracks + float d = std::pow((pp.pX - posX) * sigmaX, 2) + std::pow((pp.pZ - posZ) * sigmaZ, 2); // TODO different sigma for tracks if (d < trackdist) { trackdist = d; trackDx = pp.pX - posX; @@ -1095,10 +1100,10 @@ struct caloClusterProducerTask { } if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + cpvdist = std::sqrt(cpvdist); // was squared } if (trackdist != 99.) { // was evaluated - trackdist = sqrt(trackdist); // was squared + trackdist = std::sqrt(trackdist); // was squared } float lambdaShort = 0., lambdaLong = 0.; @@ -1145,7 +1150,7 @@ struct caloClusterProducerTask { //------------------------------------------------------------ void processFullMC(o2::aod::BCsWithTimestamps const& bcs, o2::aod::Collisions const& colls, - mcCells& cells, + McCells const& cells, o2::aod::CaloTriggers const&, o2::aod::CPVClusters const& cpvs, o2::aod::FullTracks const& tracks) @@ -1169,7 +1174,7 @@ struct caloClusterProducerTask { std::map bcMap; int bcId = 0; - for (auto bc : bcs) { + for (const auto & bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -1177,7 +1182,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (auto cl : colls) { + for (const auto & cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -1194,11 +1199,11 @@ struct caloClusterProducerTask { // Fill output table // calibration may be updated by CCDB fetcher - const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(mBadMapPath, timestamp); - const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(mCalibPath, timestamp); + const o2::phos::BadChannelsMap* badMap = ccdb->getForTimeStamp(badMapPath, timestamp); + const o2::phos::CalibParams* calibParams = ccdb->getForTimeStamp(calibPath, timestamp); if (!isMC && !skipL1phase) { - const std::vector* vec = ccdb->getForTimeStamp>(mL1PhasePath, timestamp); + const std::vector* vec = ccdb->getForTimeStamp>(l1phasePath, timestamp); if (vec) { clusterizerPHOS->setL1phase((*vec)[0]); } else { @@ -1228,7 +1233,7 @@ struct caloClusterProducerTask { o2::InteractionRecord ir; const int kPHOS = 0; - for (auto& c : cells) { + for (const auto& c : cells) { if (c.caloType() != kPHOS) // PHOS continue; // Fix for bug in trigger digits @@ -1283,7 +1288,7 @@ struct caloClusterProducerTask { // Find CPV clusters corresponding to PHOS trigger records std::vector> cpvMatchPoints[kCpvCells]; // Number of entries in each cell per TrigRecord - std::vector cpvNMatchPoints; + std::vector cpvNMatchPoints; cpvNMatchPoints.reserve(outputPHOSClusterTrigRecs.size()); int64_t curBC = -1; @@ -1312,7 +1317,7 @@ struct caloClusterProducerTask { if (cpvclu.amplitude() < static_cast>(cpvMinE)[static_cast(cpvclu.moduleNumber()) - 2]) { continue; } - int index = CpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); + int index = cpvMatchIndex(cpvclu.moduleNumber(), cpvclu.posX(), cpvclu.posZ()); cpvMatchPoints[index].emplace_back(cpvclu.posX(), cpvclu.posZ()); } if (cpvNMatchPoints.size()) { @@ -1321,9 +1326,9 @@ struct caloClusterProducerTask { } } // same for tracks - std::vector trackMatchPoints[kCpvCells]; // tracks hit in grid/cell in PHOS + std::vector trackMatchPoints[kCpvCells]; // tracks hit in grid/cell in PHOS // Number of entries in each cell per TrigRecord - std::vector trackNMatchPoints; + std::vector trackNMatchPoints; trackNMatchPoints.reserve(outputPHOSClusterTrigRecs.size()); curBC = -1; @@ -1334,7 +1339,7 @@ struct caloClusterProducerTask { } } bool keepBC = false; - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { if (cluTR.getBCData().toLong() == curBC) { keepBC = true; break; @@ -1361,7 +1366,7 @@ struct caloClusterProducerTask { curBC = track.collision().bc_as().globalBC(); } keepBC = false; - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { if (cluTR.getBCData().toLong() == curBC) { keepBC = true; break; @@ -1385,7 +1390,7 @@ struct caloClusterProducerTask { float trackX, trackZ; auto trackPar = getTrackPar(track); if (impactOnPHOS(trackPar, track.trackEtaEmcal(), track.trackPhiEmcal(), track.collision().posZ(), module, trackX, trackZ)) { - int index = CpvMatchIndex(module, trackX, trackZ); + int index = cpvMatchIndex(module, trackX, trackZ); trackMatchPoints[index].emplace_back(trackX, trackZ, track.globalIndex()); } } @@ -1396,7 +1401,7 @@ struct caloClusterProducerTask { } // Fill output tables - for (auto& cluTR : outputPHOSClusterTrigRecs) { + for (const auto& cluTR : outputPHOSClusterTrigRecs) { int firstClusterInEvent = cluTR.getFirstEntry(); int lastClusterInEvent = firstClusterInEvent + cluTR.getNumberOfObjects(); @@ -1442,7 +1447,7 @@ struct caloClusterProducerTask { // Correction for the depth of the shower starting point (TDR p 127) const float para = 0.925; const float parb = 6.52; - float depth = para * TMath::Log(e) + parb; + float depth = para * std::log(e) + parb; posX -= posX * depth / 460.; posZ -= (posZ - vtx.Z()) * depth / 460.; @@ -1455,52 +1460,52 @@ struct caloClusterProducerTask { continue; } - e = Nonlinearity(e); + e = nonlinearity(e); mom.SetMag(e); // CPV and track match - const float cellSizeX = 2 * cpvMaxX / kCpvX; - const float cellSizeZ = 2 * cpvMaxZ / kCpvZ; + const float cellSizeX = 2 * kCpvMaxX / kCpvX; + const float cellSizeZ = 2 * kCpvMaxZ / kCpvZ; // look 9 CPV regions around PHOS cluster - int phosIndex = CpvMatchIndex(mod, posX, posZ); + int phosIndex = cpvMatchIndex(mod, posX, posZ); std::vector regions; regions.push_back(phosIndex); - if (posX > -cpvMaxX + cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom left + if (posX > -kCpvMaxX + cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom left regions.push_back(phosIndex - kCpvZ - 1); } regions.push_back(phosIndex - kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top left + if (posZ < kCpvMaxZ - cellSizeZ) { // top left regions.push_back(phosIndex - kCpvZ + 1); } } - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom regions.push_back(phosIndex - 1); } - if (posZ < cpvMaxZ - cellSizeZ) { // top + if (posZ < kCpvMaxZ - cellSizeZ) { // top regions.push_back(phosIndex + 1); } - if (posX < cpvMaxX - cellSizeX) { - if (posZ > -cpvMaxZ + cellSizeZ) { // bottom right + if (posX < kCpvMaxX - cellSizeX) { + if (posZ > -kCpvMaxZ + cellSizeZ) { // bottom right regions.push_back(phosIndex + kCpvZ - 1); } regions.push_back(phosIndex + kCpvZ); - if (posZ < cpvMaxZ - cellSizeZ) { // top right + if (posZ < kCpvMaxZ - cellSizeZ) { // top right regions.push_back(phosIndex + kCpvZ + 1); } } - float sigmaX = 1. / TMath::Min(5.2, 1.111 + 0.56 * TMath::Exp(-0.031 * e * e) + 4.8 / TMath::Power(e + 0.61, 3)); // inverse sigma X - float sigmaZ = 1. / TMath::Min(3.3, 1.12 + 0.35 * TMath::Exp(-0.032 * e * e) + 0.75 / TMath::Power(e + 0.24, 3)); // inverse sigma Z + float sigmaX = 1. / std::min(5.2, 1.111 + 0.56 * std::exp(-0.031 * e * e) + 4.8 / std::pow(e + 0.61, 3)); // inverse sigma X + float sigmaZ = 1. / std::min(3.3, 1.12 + 0.35 * std::exp(-0.032 * e * e) + 0.75 / std::pow(e + 0.24, 3)); // inverse sigma Z float cpvdist = 99., trackdist = 99.; // float cpvDx = 0., cpvDz = 0.; float trackDx = 9999., trackDz = 9999.; int trackindex = -1; - for (int indx : regions) { + for (const int indx : regions) { if (cpvPoints != cpvNMatchPoints.end()) { if (indx >= 0 && indx < kCpvCells) { for (int ii = cpvPoints->mStart[indx]; ii < cpvPoints->mEnd[indx]; ii++) { auto p = cpvMatchPoints[indx][ii]; - float d = pow((p.first - posX) * sigmaX, 2) + pow((p.second - posZ) * sigmaZ, 2); + float d = std::pow((p.first - posX) * sigmaX, 2) + std::pow((p.second - posZ) * sigmaZ, 2); if (d < cpvdist) { cpvdist = d; } @@ -1511,7 +1516,7 @@ struct caloClusterProducerTask { if (trackPoints != trackNMatchPoints.end()) { for (int ii = trackPoints->mStart[indx]; ii < trackPoints->mEnd[indx]; ii++) { auto pp = trackMatchPoints[indx][ii]; - float d = pow((pp.pX - posX) * sigmaX, 2) + pow((pp.pZ - posZ) * sigmaZ, 2); // TODO different sigma for tracks + float d = std::pow((pp.pX - posX) * sigmaX, 2) + std::pow((pp.pZ - posZ) * sigmaZ, 2); // TODO different sigma for tracks if (d < trackdist) { trackdist = d; trackDx = pp.pX - posX; @@ -1523,10 +1528,10 @@ struct caloClusterProducerTask { } if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + cpvdist = std::sqrt(cpvdist); // was squared } if (trackdist != 99.) { // was evaluated - trackdist = sqrt(trackdist); // was squared + trackdist = std::sqrt(trackdist); // was squared } float lambdaShort = 0., lambdaLong = 0.; @@ -1540,7 +1545,7 @@ struct caloClusterProducerTask { mclabels.clear(); mcamplitudes.clear(); gsl::span spDigList = outputTruthCont.getLabels(i); - for (auto cellLab : spDigList) { + for (const auto & cellLab : spDigList) { mclabels.push_back(cellLab.getTrackID()); // Track ID in current event? mcamplitudes.push_back(cellLab.getEdep()); } @@ -1583,15 +1588,15 @@ struct caloClusterProducerTask { PROCESS_SWITCH(caloClusterProducerTask, processFullMC, "Process MC with track matching", false); - int CpvMatchIndex(int16_t module, float x, float z) + int cpvMatchIndex(int16_t module, float x, float z) { // calculate cell index in grid over PHOS detector - const float cellSizeX = 2 * cpvMaxX / kCpvX; - const float cellSizeZ = 2 * cpvMaxZ / kCpvZ; + const float cellSizeX = 2 * kCpvMaxX / kCpvX; + const float cellSizeZ = 2 * kCpvMaxZ / kCpvZ; // in track matching tracks can be beyond CPV surface // assign these tracks to the closest cell - int ix = std::max(0, static_cast((x + cpvMaxX) / cellSizeX)); - int iz = std::max(0, static_cast((z + cpvMaxZ) / cellSizeZ)); + int ix = std::max(0, static_cast((x + kCpvMaxX) / cellSizeX)); + int iz = std::max(0, static_cast((z + kCpvMaxZ) / cellSizeZ)); if (ix >= kCpvX) { ix = kCpvX - 1; } @@ -1612,17 +1617,12 @@ struct caloClusterProducerTask { const float etaMax = 0.178266; double bz = o2::base::Propagator::Instance()->getNominalBz(); // magnetic field - if (trackPhi < phiMin || trackPhi > phiMax || abs(trackEta) > etaMax) { // do not match even approximately + trackPhi = RecoDecay::constrainAngle(trackPhi,0.,1); //constrain angle to range 0,twoPi + if (trackPhi < phiMin || trackPhi > phiMax || std::abs(trackEta) > etaMax) { // do not match even approximately return false; } const float dphi = 20. * 0.017453293; - if (trackPhi < 0.) { - trackPhi += TMath::TwoPi(); - } - if (trackPhi > TMath::TwoPi()) { - trackPhi -= TMath::TwoPi(); - } module = 1 + static_cast((trackPhi - phiMin) / dphi); if (module < 1) { module = 1; @@ -1632,11 +1632,11 @@ struct caloClusterProducerTask { } // get PHOS radius - constexpr float shiftY = -1.26; // Depth-optimized + const double shiftY = -1.26; // Depth-optimized double posL[3] = {0., 0., shiftY}; // local position at the center of module double posG[3] = {0}; geomPHOS->getAlignmentMatrix(module)->LocalToMaster(posL, posG); - double rPHOS = sqrt(posG[0] * posG[0] + posG[1] * posG[1]); + double rPHOS = std::sqrt(posG[0] * posG[0] + posG[1] * posG[1]); double alpha = (230. + 20. * module) * 0.017453293; // During main reconstruction track was propagated to radius 460 cm with accounting material @@ -1661,7 +1661,7 @@ struct caloClusterProducerTask { return false; } alpha = trackPar.getAlpha(); - double ca = cos(alpha), sa = sin(alpha); + double ca = std::cos(alpha), sa = std::sin(alpha); posG[0] = trackPar.getX() * ca - trackPar.getY() * sa; posG[1] = trackPar.getY() * ca + trackPar.getX() * sa; posG[2] = trackPar.getZ(); @@ -1670,7 +1670,7 @@ struct caloClusterProducerTask { trackX = posL[0]; trackZ = posL[1]; // If trackX beyond the module, switch to the next one - if (abs(trackX) < xmax || (trackX < -xmax && module == 1) || (trackX > xmax && module == 4)) { + if (std::abs(trackX) < xmax || (trackX < -xmax && module == 1) || (trackX > xmax && module == 4)) { return true; } // re-do extrapolation to correct module @@ -1687,7 +1687,7 @@ struct caloClusterProducerTask { posL[1] = 0.; posL[2] = shiftY; // local position at the center of module geomPHOS->getAlignmentMatrix(module)->LocalToMaster(posL, posG); - rPHOS = sqrt(posG[0] * posG[0] + posG[1] * posG[1]); + rPHOS = std::sqrt(posG[0] * posG[0] + posG[1] * posG[1]); if (!trackPar.rotate(alpha) || !prop->PropagateToXBxByBz(trackPar, xtrg, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrNONE)) { @@ -1703,8 +1703,8 @@ struct caloClusterProducerTask { return false; } alpha = trackPar.getAlpha(); - ca = cos(alpha); - sa = sin(alpha); + ca = std::cos(alpha); + sa = std::sin(alpha); posG[0] = trackPar.getX() * ca - trackPar.getY() * sa; posG[1] = trackPar.getY() * ca + trackPar.getX() * sa; posG[2] = trackPar.getZ(); @@ -1715,13 +1715,37 @@ struct caloClusterProducerTask { return true; } - float Nonlinearity(float en) + float nonlinearity(float en) { // Correct for non-linearity - switch (mNonlinType) { + switch (nonlinType) { case 0: return en; - case 1: { + case 1: { // Data Run3 + const double a = 9.3494e-01; + const double b = 1.00526e-02; + const double c = 8.45164e-02; + const double d = -1.03364e-02; + const double f = 5.4803e-03; + const double g = 0.779983; + const double h = 0.622282; + const double k = 8.0182e-05; + double eMin = std::max(float(0.1), en); // Parameterization valid down to 100 MeV + return en * (a + b * eMin + c / eMin + d / (eMin * eMin) + f / ((eMin - g) * (eMin - g) + h * h) + k / std::pow(eMin, 4)); + } + case 2: { // MC + const double a = 1.14875; + const double b = -1.24286e-04; + const double c = -0.0498217; + const double d = -0.00215362; + const double f = 0.886539; + const double g = -1.98282; + const double h = 0.0178562; + const double k = 5.03164e-04; + double eMin = std::max(float(0.1), en); // Parameterization valid down to 100 MeV + return en * (a + b * eMin + c / eMin + d / (eMin * eMin) + f / ((eMin - g) * (eMin - g) + h * h) + k / std::pow(eMin, 4)); + } + case 3: { // Obsolete data Run3 const double a = 9.34913e-01; const double b = 2.33e-03; const double c = -8.10e-05; diff --git a/PWGEM/Tasks/CMakeLists.txt b/PWGEM/Tasks/CMakeLists.txt index 4f0a92bd6f2..e1bb6a8e85c 100644 --- a/PWGEM/Tasks/CMakeLists.txt +++ b/PWGEM/Tasks/CMakeLists.txt @@ -26,7 +26,7 @@ o2physics_add_dpl_workflow(phostrigqa o2physics_add_dpl_workflow(phospi0 SOURCES phosPi0.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(phoscalib @@ -46,7 +46,7 @@ o2physics_add_dpl_workflow(phosnbar o2physics_add_dpl_workflow(phosnonlin SOURCES phosNonlin.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase O2::PHOSReconstruction + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase O2::PHOSReconstruction O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(phoselid diff --git a/PWGEM/Tasks/phosNonlin.cxx b/PWGEM/Tasks/phosNonlin.cxx index 9115fa432c8..a0d28285326 100644 --- a/PWGEM/Tasks/phosNonlin.cxx +++ b/PWGEM/Tasks/phosNonlin.cxx @@ -9,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file phosNonlin.cxx +/// \brief task to calculate PHOS non-lienarity based on pi0 peak position +/// \author Dmitri Peresunko +/// + #include #include #include @@ -20,6 +25,8 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" #include "Framework/ConfigParamSpec.h" #include "Framework/runDataProcessing.h" @@ -34,27 +41,23 @@ #include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPLHCIFData.h" -/// \struct PHOS pi0 analysis -/// \brief Monitoring task for PHOS related quantities -/// \author Dmitri Peresunko, NRC "Kurchatov institute" -/// \since Nov, 2022 -/// - using namespace o2; using namespace o2::aod::evsel; using namespace o2::framework; using namespace o2::framework::expressions; struct phosNonlin { - Configurable mEvSelTrig{"mEvSelTrig", kTVXinPHOS, "Select events with this trigger"}; - Configurable mParamType{"mParamType", 0, "Functional form 0: a-la data, 1: a-la MC"}; - Configurable mMinCluE{"mMinCluE", 0.1, "Minimum cluster energy for analysis"}; - Configurable mMinCluTime{"minCluTime", -25.e-9, "Min. cluster time"}; - Configurable mMaxCluTime{"maxCluTime", 25.e-9, "Max. cluster time"}; - Configurable mMinCluNcell{"minCluNcell", 1, "min cells in cluster"}; - Configurable mMinM02{"minM02", 0.2, "Min disp M02 cut"}; - Configurable mNMixedEvents{"nMixedEvents", 2, "number of events to mix"}; - Configurable mSelectOneCollPerBC{"selectOneColPerBC", true, "skip multiple coll. per bc"}; + Configurable skimmedProcessing{"skimmedProcessing", true, "Skimmed dataset processing"}; + Configurable mTrigName{"trigName", "fPHOSPhoton", "name of offline trigger"}; + Configurable zorroCCDBpath{"zorroCCDBpath", "/Users/m/mpuccio/EventFiltering/OTS/", "path to the zorro ccdb objects"}; + Configurable evSelTrig{"evSelTrig", kTVXinPHOS, "Select events with this trigger"}; + Configurable paramType{"paramType", 0, "Functional form 0: a-la data, 1: a-la MC"}; + Configurable minCluE{"minCluE", 0.1, "Minimum cluster energy for analysis"}; + Configurable minCluTime{"minCluTime", -25.e-9, "Min. cluster time"}; + Configurable maxCluTime{"maxCluTime", 25.e-9, "Max. cluster time"}; + Configurable minCluNcell{"minCluNcell", 1, "min cells in cluster"}; + Configurable minM02{"minM02", 0.2, "Min disp M02 cut"}; + Configurable nMixedEvents{"nMixedEvents", 2, "number of events to mix"}; Configurable mA{"mA", 9.34913e-01, "A"}; Configurable mdAi{"mdAi", 0., "A var. vs i"}; Configurable mdAj{"mdAj", 0., "A var. vs j"}; @@ -82,25 +85,30 @@ struct phosNonlin { using SelCollisions = soa::Join; using BCsWithBcSels = soa::Join; + o2::framework::Service ccdb; HistogramRegistry mHistManager1{"phosNonlinHistograms"}; + Zorro zorro; + OutputObj zorroSummary{"zorroSummary"}; + // class to keep photon candidate parameters - class photon : public TLorentzVector + class Photon : public TLorentzVector { public: - photon() = default; - photon(const photon& p) = default; - photon(double px, double py, double pz, double e, int m) : TLorentzVector(px, py, pz, e), mod(m) {} + Photon() = default; + Photon(const Photon& p) = default; + Photon(double px, double py, double pz, double e, int m) : TLorentzVector(px, py, pz, e), mod(m) {} public: int mod; }; + int mRunNumber = -1; // Current run number int mixedEventBin = 0; // Which list of Mixed use for mixing - std::vector mCurEvent; - static constexpr int nMaxMixBins = 10; // maximal number of kinds of events for mixing - std::array>, nMaxMixBins> mMixedEvents; + std::vector mCurEvent; + static constexpr int kMaxMixBins = 10; // maximal number of kinds of events for mixing + std::array>, kMaxMixBins> mMixedEvents; // fast access to histos static constexpr int mNp = 10; @@ -116,6 +124,12 @@ struct phosNonlin { { LOG(info) << "Initializing PHOS nonlin analysis task ..."; + zorroSummary.setObject(zorro.getZorroSummary()); + zorro.setBaseCCDBPath(zorroCCDBpath.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + const AxisSpec ptAxis{pt, "p_{T} (GeV/c)"}, mggAxis{150, 0., 0.3, "m_{#gamma#gamma} (GeV/c^{2})"}; @@ -141,27 +155,41 @@ struct phosNonlin { /// \brief Process PHOS data void process(SelCollisions::iterator const& col, - aod::CaloClusters const& clusters) + aod::CaloClusters const& clusters, + aod::BCsWithTimestamps const&) { // Fill number of events of different kind - if (!col.alias_bit(mEvSelTrig)) { - return; + if (skimmedProcessing) { + auto bc = col.template bc_as(); + if (mRunNumber != bc.runNumber()) { + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), trigName); + zorro.populateHistRegistry(mHistManager1, bc.runNumber()); + mRunNumber = bc.runNumber(); + } + + if (!zorro.isSelected(bc.globalBC())) { + return; /// + } + } else { + if (!col.selection_bit(evSelTrig)) { + return; + } } mixedEventBin = findMixedEventBin(col.posZ()); mCurEvent.clear(); int i, j, k, l; - for (const auto& clu : clusters) { - if (clu.e() < mMinCluE || - clu.ncell() < mMinCluNcell || - clu.time() > mMaxCluTime || clu.time() < mMinCluTime || - clu.m02() < mMinM02) { + for (auto const & clu : clusters) { + if (clu.e() < minCluE || + clu.ncell() < minCluNcell || + clu.time() > maxCluTime || clu.time() < minCluTime || + clu.m02() < minM02) { continue; } - photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod()); + Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod()); // Mix with other photons added to stack - for (auto ph2 : mCurEvent) { + for (auto const ph2 : mCurEvent) { double m = (ph1 + ph2).M(); double pt1 = ph1.Pt(); double pt2 = ph2.Pt(); @@ -170,8 +198,8 @@ struct phosNonlin { l = mNp / 2; for (i = 0; i < mNp; i++) { for (j = 0; j < mNp; j++) { - if (ph1.E() * NonLin(ph1.E(), i, j, k, l) > mMinCluE && ph2.E() * NonLin(ph2.E(), i, j, k, l) > mMinCluE) { - Double_t m12 = m * TMath::Sqrt(NonLin(ph1.E(), i, j, k, l) * NonLin(ph2.E(), i, j, k, l)); + if (ph1.E() * nonLin(ph1.E(), i, j, k, l) > minCluE && ph2.E() * nonLin(ph2.E(), i, j, k, l) > minCluE) { + double m12 = m * std::sqrt(nonLin(ph1.E(), i, j, k, l) * nonLin(ph2.E(), i, j, k, l)); hReIJ[i * mNp + j]->Fill(m12, pt1); hReIJ[i * mNp + j]->Fill(m12, pt2); // if(ph1.mod==ph2.mod){ @@ -185,8 +213,8 @@ struct phosNonlin { j = mNp / 2; for (k = 0; k < mNp; k++) { for (l = 0; l < mNp; l++) { - if (ph1.E() * NonLin(ph1.E(), i, j, k, l) > mMinCluE && ph2.E() * NonLin(ph2.E(), i, j, k, l) > mMinCluE) { - Double_t m12 = m * TMath::Sqrt(NonLin(ph1.E(), i, j, k, l) * NonLin(ph2.E(), i, j, k, l)); + if (ph1.E() * nonLin(ph1.E(), i, j, k, l) > minCluE && ph2.E() * nonLin(ph2.E(), i, j, k, l) > minCluE) { + double m12 = m * std::sqrt(nonLin(ph1.E(), i, j, k, l) * nonLin(ph2.E(), i, j, k, l)); hReKL[k * mNp + l]->Fill(m12, pt1); hReKL[k * mNp + l]->Fill(m12, pt2); // if(ph1.mod==ph2.mod){ @@ -202,9 +230,9 @@ struct phosNonlin { } // Mixed - for (auto ph1 : mCurEvent) { - for (auto mixEvent : mMixedEvents[mixedEventBin]) { - for (auto ph2 : mixEvent) { + for (auto const ph1 : mCurEvent) { + for (auto const mixEvent : mMixedEvents[mixedEventBin]) { + for (auto const ph2 : mixEvent) { double m = (ph1 + ph2).M(); double pt1 = ph1.Pt(); double pt2 = ph2.Pt(); @@ -213,8 +241,8 @@ struct phosNonlin { j = mNp / 2; k = mNp / 2; l = mNp / 2; - if (ph1.E() * NonLin(ph1.E(), i, j, k, l) > mMinCluE && ph2.E() * NonLin(ph2.E(), i, j, k, l) > mMinCluE) { - Double_t m12 = m * TMath::Sqrt(NonLin(ph1.E(), i, j, k, l) * NonLin(ph2.E(), i, j, k, l)); + if (ph1.E() * nonLin(ph1.E(), i, j, k, l) > minCluE && ph2.E() * nonLin(ph2.E(), i, j, k, l) > minCluE) { + double m12 = m * std::sqrt(nonLin(ph1.E(), i, j, k, l) * nonLin(ph2.E(), i, j, k, l)); hMi->Fill(m12, pt1); hMi->Fill(m12, pt2); } @@ -225,7 +253,7 @@ struct phosNonlin { // Fill events to store and remove oldest to keep buffer size if (mCurEvent.size() > 0) { mMixedEvents[mixedEventBin].emplace_back(mCurEvent); - if (mMixedEvents[mixedEventBin].size() > static_cast(mNMixedEvents)) { + if (mMixedEvents[mixedEventBin].size() > static_cast(nMixedEvents)) { mMixedEvents[mixedEventBin].pop_front(); } } @@ -240,33 +268,33 @@ struct phosNonlin { if (res < 0) return 0; - if (res >= nMaxMixBins) - return nMaxMixBins - 1; + if (res >= kMaxMixBins) + return kMaxMixBins - 1; return res; } //_____________________________________________________________________________ - double NonLin(double en, int i, int j, int k, int l) + double nonLin(double en, int i, int j, int k, int l) { if (en <= 0.) return 0.; - if (mParamType == 0) { - const Double_t a = mA + mdAi * (i - mNp / 2) + mdAj * (j - mNp / 2); - const Double_t b = mB + mdBi * (i - mNp / 2) + mdBj * (j - mNp / 2); - const Double_t c = mC + mdCi * (i - mNp / 2) + mdCj * (j - mNp / 2); - const Double_t d = mD + mdDk * (k - mNp / 2) + mdDl * (l - mNp / 2); - const Double_t e = mE + mdEk * (k - mNp / 2) + mdEl * (l - mNp / 2); - const Double_t f = mF + mdFk * (k - mNp / 2) + mdFl * (l - mNp / 2); - const Double_t g = mG + mdGk * (k - mNp / 2) + mdGl * (l - mNp / 2); - const Double_t s = mS + mdSi * (i - mNp / 2) + mdSj * (j - mNp / 2); + if (paramType == 0) { + const double a = mA + mdAi * (i - mNp / 2) + mdAj * (j - mNp / 2); + const double b = mB + mdBi * (i - mNp / 2) + mdBj * (j - mNp / 2); + const double c = mC + mdCi * (i - mNp / 2) + mdCj * (j - mNp / 2); + const double d = mD + mdDk * (k - mNp / 2) + mdDl * (l - mNp / 2); + const double e = mE + mdEk * (k - mNp / 2) + mdEl * (l - mNp / 2); + const double f = mF + mdFk * (k - mNp / 2) + mdFl * (l - mNp / 2); + const double g = mG + mdGk * (k - mNp / 2) + mdGl * (l - mNp / 2); + const double s = mS + mdSi * (i - mNp / 2) + mdSj * (j - mNp / 2); return a + b / en + c / (en * en) + d / ((en - e) * (en - e) + f * f) + g * en + s / (en * en * en * en); } - if (mParamType == 1) { - const Double_t a = mA + mdAi * (i - mNp / 2) + mdAj * (j - mNp / 2); - const Double_t b = mB + mdBi * (i - mNp / 2) + mdBj * (j - mNp / 2); - const Double_t c = mC + mdCi * (i - mNp / 2) + mdCj * (j - mNp / 2); - const Double_t d = mD + mdDk * (k - mNp / 2) + mdDl * (l - mNp / 2); - const Double_t e = mE + mdEk * (k - mNp / 2) + mdEl * (l - mNp / 2); - return a + b / TMath::Sqrt(en) + c / en + d / (en * TMath::Sqrt(en)) + e / (en * en); + if (paramType == 1) { + const double a = mA + mdAi * (i - mNp / 2) + mdAj * (j - mNp / 2); + const double b = mB + mdBi * (i - mNp / 2) + mdBj * (j - mNp / 2); + const double c = mC + mdCi * (i - mNp / 2) + mdCj * (j - mNp / 2); + const double d = mD + mdDk * (k - mNp / 2) + mdDl * (l - mNp / 2); + const double e = mE + mdEk * (k - mNp / 2) + mdEl * (l - mNp / 2); + return a + b / std::sqrt(en) + c / en + d / (en * std::sqrt(en)) + e / (en * en); } return 0.; } diff --git a/PWGEM/Tasks/phosPi0.cxx b/PWGEM/Tasks/phosPi0.cxx index fe248dbeaa6..77aeb76f081 100644 --- a/PWGEM/Tasks/phosPi0.cxx +++ b/PWGEM/Tasks/phosPi0.cxx @@ -9,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file phosPi0.cxx +/// \brief PHOS pi0 analysis +/// \author Dmitri Peresunko +/// + #include #include #include @@ -27,53 +32,54 @@ #include "Framework/ASoAHelpers.h" #include "Framework/HistogramRegistry.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" + #include "PHOSBase/Geometry.h" #include "CommonDataFormat/InteractionRecord.h" #include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPLHCIFData.h" -/// \struct PHOS pi0 analysis -/// \brief Monitoring task for PHOS related quantities -/// \author Dmitri Peresunko, NRC "Kurchatov institute" -/// \since Nov, 2022 -/// - using namespace o2; using namespace o2::aod::evsel; using namespace o2::framework; using namespace o2::framework::expressions; struct phosPi0 { - Configurable mIsMC{"isMC", false, "to fill MC histograms"}; - Configurable mEvSelTrig{"mEvSelTrig", kTVXinPHOS, "Select events with this trigger"}; - Configurable mMinCluE{"mMinCluE", 0.3, "Minimum cluster energy for analysis"}; - Configurable mMinCluTime{"minCluTime", -25.e-9, "Min. cluster time"}; - Configurable mMaxCluTime{"maxCluTime", 25.e-9, "Max. cluster time"}; - Configurable mMinCluNcell{"minCluNcell", 2, "min cells in cluster"}; - Configurable mMinM02{"minM02", 0.2, "Min disp M02 cut"}; - Configurable mCPVCut{"CPVCut", 2., "Min distance to track"}; - Configurable mNMixedEvents{"nMixedEvents", 10, "number of events to mix"}; - Configurable mSelectOneCollPerBC{"selectOneColPerBC", true, "skip multiple coll. per bc"}; - Configurable mFillQC{"fillQC", true, "Fill QC histos"}; - Configurable mOccE{"minOccE", 0.5, "Min. cluster energy of occupancy plots"}; + Configurable skimmedProcessing{"skimmedProcessing", false, "Skimmed dataset processing"}; + Configurable trigName{"trigName", "fPHOSPhoton", "name of offline trigger"}; + Configurable zorroCCDBpath{"zorroCCDBpath", "/Users/m/mpuccio/EventFiltering/OTS/", "path to the zorro ccdb objects"}; + Configurable evSelTrig{"evSelTrig", aod::evsel::kIsTriggerTVX, "Select events with this trigger"}; + Configurable isMC{"isMC", false, "to fill MC histograms"}; + Configurable minCluE{"minCluE", 0.3, "Minimum cluster energy for analysis"}; + Configurable minCluTime{"minCluTime", -25.e-9, "Min. cluster time"}; + Configurable maxCluTime{"maxCluTime", 25.e-9, "Max. cluster time"}; + Configurable minCluNcell{"minCluNcell", 2, "min cells in cluster"}; + Configurable minM02{"minM02", 0.2, "Min disp M02 cut"}; + Configurable cpvCut{"cpvCut", 2., "Min distance to track"}; + Configurable nMixedEvents{"nMixedEvents", 10, "number of events to mix"}; + Configurable fillQC{"fillQC", true, "Fill QC histos"}; + Configurable minOccE{"minOccE", 0.5, "Min. cluster energy of occupancy plots"}; using SelCollisions = soa::Join; using SelCollisionsMC = soa::Join; using BCsWithBcSels = soa::Join; - using mcClusters = soa::Join; - using mcAmbClusters = soa::Join; + using McClusters = soa::Join; + using McAmbClusters = soa::Join; o2::framework::Service ccdb; + Zorro zorro; + OutputObj zorroSummary{"zorroSummary"}; HistogramRegistry mHistManager{"phosPi0Histograms"}; // class to keep photon candidate parameters - class photon + class Photon { public: - photon() = default; - photon(double x, double y, double z, double ee, int m, bool isDispOK, bool isCPVOK, int mcLabel) : px(x), py(y), pz(z), e(ee), mod(m), mPID(isDispOK << 1 | isCPVOK << 2), label(mcLabel) {} - ~photon() = default; + Photon() = default; + Photon(double x, double y, double z, double ee, int m, bool isDispOK, bool isCPVOK, int mcLabel) : px(x), py(y), pz(z), e(ee), mod(m), mPID(isDispOK << 1 | isCPVOK << 2), label(mcLabel) {} + ~Photon() = default; bool isCPVOK() { return (mPID >> 2) & 1; } bool isDispOK() { return (mPID >> 1) & 1; } @@ -88,11 +94,12 @@ struct phosPi0 { int label = -1; // label of MC particle }; + int mRunNumber = 0; // Current run number int mixedEventBin = 0; // Which list of Mixed use for mixing - std::vector mCurEvent; - static constexpr int nMaxMixBins = 20; // maximal number of kinds of events for mixing - std::array>, nMaxMixBins> mMixedEvents; - std::array>, nMaxMixBins> mAmbMixedEvents; + std::vector mCurEvent; + static constexpr int kMaxMixBins = 20; // maximal number of kinds of events for mixing + std::array>, kMaxMixBins> mMixedEvents; + std::array>, kMaxMixBins> mAmbMixedEvents; int mPrevMCColId = -1; // mark MC collissions already scanned // fast access to histos @@ -110,13 +117,15 @@ struct phosPi0 { void init(InitContext const&) { LOG(info) << "Initializing PHOS pi0 analysis task ..."; + zorroSummary.setObject(zorro.getZorroSummary()); + zorro.setBaseCCDBPath(zorroCCDBpath.value); const AxisSpec ptAxis{pt, "p_{T} (GeV/c)"}, mggAxis{625, 0., 1.25, "m_{#gamma#gamma} (GeV/c^{2})"}, timeAxis{100, -500.e-9, 500.e-9, "t (s)"}, - M02Axis{100, 0., 20., "M02 (cm^{2})"}, - M20Axis{100, 0., 20., "M20 (cm^{2})"}, + m02Axis{100, 0., 20., "M02 (cm^{2})"}, + m20Axis{100, 0., 20., "M20 (cm^{2})"}, nCellsAxis{100, 0., 100., "N_{cell}"}, zAxis{56, -63., 63., "z", "z (cm)"}, phiAxis{64, -72., 72., "x", "x (cm)"}, @@ -127,15 +136,16 @@ struct phosPi0 { centAxis{10, 0., 10.}, centralityAxis{100, 0., 100., "centrality", "centrality"}; - hColl = std::get>(mHistManager.add("eventsCol", "Number of events", HistType::kTH1F, {{9, 0., 9.}})).get(); + hColl = std::get>(mHistManager.add("eventsCol", "Number of events", HistType::kTH1F, {{10, 0., 10.}})).get(); hColl->GetXaxis()->SetBinLabel(1, "All"); - hColl->GetXaxis()->SetBinLabel(2, "T0a||T0c"); - hColl->GetXaxis()->SetBinLabel(3, "T0a&&T0c"); - hColl->GetXaxis()->SetBinLabel(4, "kTVXinPHOS"); - hColl->GetXaxis()->SetBinLabel(5, "kIsTriggerTVX"); - hColl->GetXaxis()->SetBinLabel(6, "PHOSClu"); - hColl->GetXaxis()->SetBinLabel(7, "PHOSClu&&kTVXinPHOS"); - hColl->GetXaxis()->SetBinLabel(8, "Accepted"); + hColl->GetXaxis()->SetBinLabel(2, "SwTr"); + hColl->GetXaxis()->SetBinLabel(3, "T0a||T0c"); + hColl->GetXaxis()->SetBinLabel(4, "T0a&&T0c"); + hColl->GetXaxis()->SetBinLabel(5, "kTVXinPHOS"); + hColl->GetXaxis()->SetBinLabel(6, "kIsTriggerTVX"); + hColl->GetXaxis()->SetBinLabel(7, "PHOSClu"); + hColl->GetXaxis()->SetBinLabel(8, "PHOSClu&&kTVXinPHOS"); + hColl->GetXaxis()->SetBinLabel(9, "Accepted"); auto h2{std::get>(mHistManager.add("eventsBC", "Number of events per trigger", HistType::kTH1F, {{8, 0., 8.}}))}; h2->GetXaxis()->SetBinLabel(1, "All"); @@ -150,14 +160,14 @@ struct phosPi0 { mHistManager.add("contributors", "Contributors per collision", HistType::kTH2F, {{10, 0., 10.}, {10, 0., 100.}}); mHistManager.add("vertex", "vertex", HistType::kTH1F, {vertexAxis}); - if (mFillQC) { + if (fillQC) { // QC histograms for normal collisions mHistManager.add("cluSp", "Cluster spectrum per module", HistType::kTH2F, {ptAxis, modAxis}); mHistManager.add("cluSpDisp", "Cluster spectrum per module", HistType::kTH2F, {ptAxis, modAxis}); mHistManager.add("cluSpCPV", "Cluster spectrum per module", HistType::kTH2F, {ptAxis, modAxis}); mHistManager.add("cluSpBoth", "Cluster spectrum per module", HistType::kTH2F, {ptAxis, modAxis}); - mHistManager.add("hM02Clu", "(M02,M20) in clusters", HistType::kTH2F, {ptAxis, M02Axis}); - mHistManager.add("hM20Clu", "(M02,M20) in clusters", HistType::kTH2F, {ptAxis, M20Axis}); + mHistManager.add("hM02Clu", "(M02,M20) in clusters", HistType::kTH2F, {ptAxis, m02Axis}); + mHistManager.add("hM20Clu", "(M02,M20) in clusters", HistType::kTH2F, {ptAxis, m20Axis}); mHistManager.add("hNcellClu", "Number of cells in clusters", HistType::kTH3F, {ptAxis, nCellsAxis, modAxis}); mHistManager.add("cluETime", "Cluster time vs E", HistType::kTH3F, {ptAxis, timeAxis, modAxis}); @@ -183,7 +193,7 @@ struct phosPi0 { HistType::kTH2F, {mggAxis, ptAxis})) .get(); - if (mIsMC) { + if (isMC) { hSignalAll = std::get>(mHistManager.add("mggSignal", "inv mass for correlated pairs", HistType::kTH2F, {mggAxis, ptAxis})) .get(); @@ -216,27 +226,49 @@ struct phosPi0 { hMiBoth = std::get>(mHistManager.add("mggMiBoth", "inv mass for centrality", HistType::kTH2F, {mggAxis, ptAxis})) .get(); - if (mIsMC) { + if (isMC) { mHistManager.add("hMCPi0SpAll", "pi0 spectrum inclusive", HistType::kTH1F, {ptAxis}); mHistManager.add("hMCPi0SpPrim", "pi0 spectrum Primary", HistType::kTH1F, {ptAxis}); mHistManager.add("hMCPi0RapPrim", "pi0 rapidity primary", HistType::kTH1F, {{100, -1., 1., "Rapidity"}}); - mHistManager.add("hMCPi0PhiPrim", "pi0 phi primary", HistType::kTH1F, {{100, 0., TMath::TwoPi(), "#phi (rad)"}}); - mHistManager.add("hMCPi0SecVtx", "pi0 secondary", HistType::kTH2F, {{100, 0., 500., "R (cm)"}, {100, -TMath::Pi(), TMath::Pi(), "#phi (rad)"}}); + mHistManager.add("hMCPi0PhiPrim", "pi0 phi primary", HistType::kTH1F, {{100, 0., o2::constants::TwoPI, "#phi (rad)"}}); + mHistManager.add("hMCPi0SecVtx", "pi0 secondary", HistType::kTH2F, {{100, 0., 500., "R (cm)"}, {100, -o2::constants::PI, o2::constants::PI, "#phi (rad)"}}); } } + // template + // float getCentrality(Tcoll const& collision) + // { + // float centrality = 1.; + // if constexpr (std::is_same::value || std::is_same::value || std::is_same::value) { + // if (cfgCentralityEstimator == nuclei::centDetectors::kFV0A) { + // centrality = collision.centFV0A(); + // } else if (cfgCentralityEstimator == nuclei::centDetectors::kFT0M) { + // centrality = collision.centFT0M(); + // } else if (cfgCentralityEstimator == nuclei::centDetectors::kFT0A) { + // centrality = collision.centFT0A(); + // } else if (cfgCentralityEstimator == nuclei::centDetectors::kFT0C) { + // centrality = collision.centFT0C(); + // } else { + // LOG(warning) << "Centrality estimator not valid. Possible values: (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3). Centrality set to 1."; + // } + // } + // return centrality; + // } + /// \brief Process PHOS data void processData(SelCollisions::iterator const& col, - aod::CaloClusters const& clusters) + aod::CaloClusters const& clusters, + aod::BCsWithTimestamps const&) { aod::McParticles const* mcPart = nullptr; scanAll(col, clusters, mcPart); } PROCESS_SWITCH(phosPi0, processData, "processData", true); void processMC(SelCollisionsMC::iterator const& col, - mcClusters const& clusters, + McClusters const& clusters, aod::McParticles const& mcPart, - aod::McCollisions const& /*mcCol*/) + aod::McCollisions const& /*mcCol*/, + aod::BCsWithTimestamps const&) { scanAll(col, clusters, &mcPart); } @@ -247,43 +279,94 @@ struct phosPi0 { TClusters& clusters, aod::McParticles const* mcPart) { - bool isColSelected = false; mixedEventBin = 0; - hColl->Fill(0.5); + if (skimmedProcessing) { + auto bc = col.template bc_as(); + if (mRunNumber != bc.runNumber()) { + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), trigName); + zorro.populateHistRegistry(mHistManager, bc.runNumber()); + mRunNumber = bc.runNumber(); + } + + if (!zorro.isSelected(bc.globalBC())) { + return; /// + } + } else { + if (!col.selection_bit(evSelTrig)) { + return; + } + } + + hColl->Fill(1.5); + const double vtxCut = 10.; + double vtxZ = col.posZ(); + mHistManager.fill(HIST("vertex"), vtxZ); + bool isColSelected = false; + if constexpr (isMC) { + isColSelected = (col.selection_bit(kIsTriggerTVX) && (clusters.size() > 0)); + } else { + isColSelected = col.selection_bit(evSelTrig) && sta::abs(vtxZ) < vtxCut; // col.alias_bit(evSelTrig) + // collision.selection_bit(aod::evsel::kNoTimeFrameBorder); + } + if (col.selection_bit(kIsBBT0A) || col.selection_bit(kIsBBT0C)) { - hColl->Fill(1.5); + hColl->Fill(2.5); } if (col.selection_bit(kIsBBT0A) && col.selection_bit(kIsBBT0C)) { - hColl->Fill(2.5); + hColl->Fill(3.5); } if (col.alias_bit(kTVXinPHOS)) { - hColl->Fill(3.5); + hColl->Fill(4.5); } if (col.selection_bit(kIsTriggerTVX)) { - hColl->Fill(4.5); + hColl->Fill(5.5); } if (clusters.size() > 0) { - hColl->Fill(5.5); + hColl->Fill(6.5); if (col.alias_bit(kTVXinPHOS)) { - hColl->Fill(6.5); + hColl->Fill(7.5); } } - isColSelected = false; - if constexpr (isMC) { - isColSelected = (col.selection_bit(kIsTriggerTVX) && (clusters.size() > 0)); - } else { - isColSelected = col.alias_bit(mEvSelTrig); - } - double vtxZ = col.posZ(); - mHistManager.fill(HIST("vertex"), vtxZ); + // //Event Plane| jet orientation + // if (flag & (kProton | kDeuteron | kTriton | kHe3 | kHe4) || doprocessMC) { /// ignore PID pre-selections for the MC + // if constexpr (std::is_same::value) { + // nuclei::candidates_flow.emplace_back(NucleusCandidateFlow{ + // collision.centFV0A(), + // collision.centFT0M(), + // collision.centFT0A(), + // collision.centFT0C(), + // collision.psiFT0A(), + // collision.multFT0A(), + // collision.psiFT0C(), + // collision.multFT0C(), + // collision.psiTPC(), + // collision.psiTPCL(), + // collision.psiTPCR(), + // collision.multTPC()}); + // } else if constexpr (std::is_same::value) { + // nuclei::candidates_flow.emplace_back(NucleusCandidateFlow{ + // collision.centFV0A(), + // collision.centFT0M(), + // collision.centFT0A(), + // collision.centFT0C(), + // 0.5 * std::atan2(collision.qvecFT0AIm(), collision.qvecFT0ARe()), + // collision.multFT0A(), + // 0.5 * std::atan2(collision.qvecFT0CIm(), collision.qvecFT0CRe()), + // collision.multFT0C(), + // -999., + // 0.5 * std::atan2(collision.qvecBNegIm(), collision.qvecBNegRe()), + // 0.5 * std::atan2(collision.qvecBPosIm(), collision.qvecBPosRe()), + // collision.multTPC()}); + // } + int mult = 1.; // multiplicity TODO!!! mixedEventBin = findMixedEventBin(vtxZ, mult); if (!isColSelected) { return; } - hColl->Fill(7.5); + hColl->Fill(8.5); // Fill MC distributions // pion rapidity, pt, phi @@ -291,7 +374,7 @@ struct phosPi0 { if constexpr (isMC) { // check current collision Id for clusters int cluMcBCId = -1; - for (auto clu : clusters) { + for (const auto &clu : clusters) { auto mcList = clu.labels(); // const std::vector int nParents = mcList.size(); for (int iParent = 0; iParent < nParents; iParent++) { // Not found nbar parent yiet @@ -309,19 +392,19 @@ struct phosPi0 { if (mcPart->begin() != mcPart->end()) { if (mcPart->begin().mcCollisionId() != mPrevMCColId) { mPrevMCColId = mcPart->begin().mcCollisionId(); // to avoid scanning full MC table each BC - for (auto part : *mcPart) { + for (const auto part : *mcPart) { if (part.mcCollision().bcId() != cluMcBCId) { continue; } if (part.pdgCode() == 111) { - double r = sqrt(pow(part.vx(), 2) + pow(part.vy(), 2)); + double r = std::sqrt(std::pow(part.vx(), 2) + std::pow(part.vy(), 2)); if (r < 0.5) { mHistManager.fill(HIST("hMCPi0RapPrim"), part.y()); } - if (abs(part.y()) < .5) { + if (std::abs(part.y()) < .5) { double pt = part.pt(); mHistManager.fill(HIST("hMCPi0SpAll"), pt); - double phiVtx = atan2(part.vy(), part.vx()); + double phiVtx = std::atan2(part.vy(), part.vx()); if (r > 0.5) { mHistManager.fill(HIST("hMCPi0SecVtx"), r, phiVtx); } else { @@ -339,31 +422,31 @@ struct phosPi0 { mCurEvent.clear(); for (const auto& clu : clusters) { // Fill QC histos - if (mFillQC) { + if (fillQC) { mHistManager.fill(HIST("hM02Clu"), clu.e(), clu.m02()); mHistManager.fill(HIST("hM20Clu"), clu.e(), clu.m20()); mHistManager.fill(HIST("hNcellClu"), clu.e(), clu.ncell(), clu.mod()); mHistManager.fill(HIST("cluETime"), clu.e(), clu.time(), clu.mod()); } - if (clu.e() < mMinCluE || - clu.ncell() < mMinCluNcell || - clu.time() > mMaxCluTime || clu.time() < mMinCluTime || - clu.m02() < mMinM02) { + if (clu.e() < minCluE || + clu.ncell() < minCluNcell || + clu.time() > maxCluTime || clu.time() < minCluTime || + clu.m02() < minM02) { continue; } - if (mFillQC) { + if (fillQC) { mHistManager.fill(HIST("cluSp"), clu.e(), clu.mod()); - if (clu.e() > mOccE) { + if (clu.e() > minOccE) { mHistManager.fill(HIST("cluOcc"), clu.x(), clu.z(), clu.mod()); if (clu.trackdist() > 2.) { mHistManager.fill(HIST("cluCPVOcc"), clu.x(), clu.z(), clu.mod()); mHistManager.fill(HIST("cluSpCPV"), clu.e(), clu.mod()); - if (TestLambda(clu.e(), clu.m02(), clu.m20())) { + if (testLambda(clu.e(), clu.m02(), clu.m20())) { mHistManager.fill(HIST("cluBothOcc"), clu.x(), clu.z(), clu.mod()); mHistManager.fill(HIST("cluSpBoth"), clu.e(), clu.mod()); } } - if (TestLambda(clu.e(), clu.m02(), clu.m20())) { + if (testLambda(clu.e(), clu.m02(), clu.m20())) { mHistManager.fill(HIST("cluDispOcc"), clu.x(), clu.z(), clu.mod()); mHistManager.fill(HIST("cluSpDisp"), clu.e(), clu.mod()); } @@ -377,17 +460,17 @@ struct phosPi0 { mcLabel = mcList[0]; } } - photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), TestLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > mCPVCut, mcLabel); + Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel); // Mix with other photons added to stack - for (auto ph2 : mCurEvent) { - double m = pow(ph1.e + ph2.e, 2) - pow(ph1.px + ph2.px, 2) - - pow(ph1.py + ph2.py, 2) - pow(ph1.pz + ph2.pz, 2); + for (cosnt auto & ph2 : mCurEvent) { + double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - + std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { - m = sqrt(m); + m = std::sqrt(m); } - double pt = sqrt(pow(ph1.px + ph2.px, 2) + - pow(ph1.py + ph2.py, 2)); - int modComb = ModuleCombination(ph1.mod, ph2.mod); + double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) + + std::pow(ph1.py + ph2.py, 2)); + int modComb = moduleCombination(ph1.mod, ph2.mod); hReMod->Fill(m, pt, modComb); hReAll->Fill(m, pt); bool isPi0 = false; @@ -427,17 +510,17 @@ struct phosPi0 { } // Mixed - for (auto ph1 : mCurEvent) { - for (auto mixEvent : mMixedEvents[mixedEventBin]) { - for (auto ph2 : mixEvent) { - double m = pow(ph1.e + ph2.e, 2) - pow(ph1.px + ph2.px, 2) - - pow(ph1.py + ph2.py, 2) - pow(ph1.pz + ph2.pz, 2); + for (const auto & ph1 : mCurEvent) { + for (const auto & mixEvent : mMixedEvents[mixedEventBin]) { + for (const auto &ph2 : mixEvent) { + double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - + std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { - m = sqrt(m); + m = std::sqrt(m); } - double pt = sqrt(pow(ph1.px + ph2.px, 2) + - pow(ph1.py + ph2.py, 2)); - int modComb = ModuleCombination(ph1.mod, ph2.mod); + double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) + + std::pow(ph1.py + ph2.py, 2)); + int modComb = moduleCombination(ph1.mod, ph2.mod); hMiMod->Fill(m, pt, modComb); hMiAll->Fill(m, pt); if (ph1.isCPVOK() && ph2.isCPVOK()) { @@ -456,7 +539,7 @@ struct phosPi0 { // Fill events to store and remove oldest to keep buffer size if (mCurEvent.size() > 0) { mMixedEvents[mixedEventBin].emplace_back(mCurEvent); - if (mMixedEvents[mixedEventBin].size() > static_cast(mNMixedEvents)) { + if (mMixedEvents[mixedEventBin].size() > static_cast(nMixedEvents)) { mMixedEvents[mixedEventBin].pop_front(); } } @@ -480,31 +563,31 @@ struct phosPi0 { mCurEvent.clear(); for (const auto& clu : clusters) { // Fill QC histos - if (mFillQC) { + if (fillQC) { mHistManager.fill(HIST("hM02Clu"), clu.e(), clu.m02()); mHistManager.fill(HIST("hM20Clu"), clu.e(), clu.m20()); mHistManager.fill(HIST("hNcellClu"), clu.e(), clu.ncell(), clu.mod()); mHistManager.fill(HIST("cluETime"), clu.e(), clu.time(), clu.mod()); } - if (clu.e() < mMinCluE || - clu.ncell() < mMinCluNcell || - clu.time() > mMaxCluTime || clu.time() < mMinCluTime || - clu.m02() < mMinM02) { + if (clu.e() < minCluE || + clu.ncell() < minCluNcell || + clu.time() > maxCluTime || clu.time() < minCluTime || + clu.m02() < minM02) { continue; } - if (mFillQC) { + if (fillQC) { mHistManager.fill(HIST("cluSp"), clu.e(), clu.mod()); - if (clu.e() > mOccE) { + if (clu.e() > minOccE) { mHistManager.fill(HIST("cluOcc"), clu.x(), clu.z(), clu.mod()); if (clu.trackdist() > 2.) { mHistManager.fill(HIST("cluCPVOcc"), clu.x(), clu.z(), clu.mod()); mHistManager.fill(HIST("cluSpCPV"), clu.e(), clu.mod()); - if (TestLambda(clu.e(), clu.m02(), clu.m20())) { + if (testLambda(clu.e(), clu.m02(), clu.m20())) { mHistManager.fill(HIST("cluBothOcc"), clu.x(), clu.z(), clu.mod()); mHistManager.fill(HIST("cluSpBoth"), clu.e(), clu.mod()); } } - if (TestLambda(clu.e(), clu.m02(), clu.m20())) { + if (testLambda(clu.e(), clu.m02(), clu.m20())) { mHistManager.fill(HIST("cluDispOcc"), clu.x(), clu.z(), clu.mod()); mHistManager.fill(HIST("cluSpDisp"), clu.e(), clu.mod()); } @@ -512,17 +595,17 @@ struct phosPi0 { } int mcLabel = -1; - photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), TestLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > mCPVCut, mcLabel); + Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel); // Mix with other photons added to stack - for (auto ph2 : mCurEvent) { - double m = pow(ph1.e + ph2.e, 2) - pow(ph1.px + ph2.px, 2) - - pow(ph1.py + ph2.py, 2) - pow(ph1.pz + ph2.pz, 2); + for (const auto & ph2 : mCurEvent) { + double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - + std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { - m = sqrt(m); + m = std::sqrt(m); } - double pt = sqrt(pow(ph1.px + ph2.px, 2) + - pow(ph1.py + ph2.py, 2)); - int modComb = ModuleCombination(ph1.mod, ph2.mod); + double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) + + std::pow(ph1.py + ph2.py, 2)); + int modComb = moduleCombination(ph1.mod, ph2.mod); hReMod->Fill(m, pt, modComb); hReAll->Fill(m, pt); @@ -542,17 +625,17 @@ struct phosPi0 { } // Mixed - for (auto ph1 : mCurEvent) { - for (auto mixEvent : mMixedEvents[mixedEventBin]) { - for (auto ph2 : mixEvent) { - double m = pow(ph1.e + ph2.e, 2) - pow(ph1.px + ph2.px, 2) - - pow(ph1.py + ph2.py, 2) - pow(ph1.pz + ph2.pz, 2); + for (const auto & ph1 : mCurEvent) { + for (const auto & mixEvent : mMixedEvents[mixedEventBin]) { + for (const auto &ph2 : mixEvent) { + double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - + std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { - m = sqrt(m); + m = std::sqrt(m); } - double pt = sqrt(pow(ph1.px + ph2.px, 2) + - pow(ph1.py + ph2.py, 2)); - int modComb = ModuleCombination(ph1.mod, ph2.mod); + double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) + + std::pow(ph1.py + ph2.py, 2)); + int modComb = moduleCombination(ph1.mod, ph2.mod); hMiMod->Fill(m, pt, modComb); hMiAll->Fill(m, pt); if (ph1.isCPVOK() && ph2.isCPVOK()) { @@ -571,7 +654,7 @@ struct phosPi0 { // Fill events to store and remove oldest to keep buffer size if (mCurEvent.size() > 0) { mMixedEvents[mixedEventBin].emplace_back(mCurEvent); - if (mMixedEvents[mixedEventBin].size() > static_cast(mNMixedEvents)) { + if (mMixedEvents[mixedEventBin].size() > static_cast(nMixedEvents)) { mMixedEvents[mixedEventBin].pop_front(); } } @@ -579,33 +662,33 @@ struct phosPi0 { PROCESS_SWITCH(phosPi0, processBC, "processBC", false); //_____________________________________________________________________________ - int ModuleCombination(int m1, int m2) + int moduleCombination(int m1, int m2) { // enumerates possible module combinations // (1,1)=0, (2,2)=1, (3,3)=2, (4,4)=3, (1,2)=(2,1)=4, (2,3)=(3,2)=5, (3,4)=(4,3)=6, (1,3)=(3,1)=7, // (2,4)=(4,2)=8, (1,4)=(4,1)=9 - int d = TMath::Abs(m1 - m2); + int d = std::abs(m1 - m2); if (d == 0) { return m1 - 1; } if (d == 1) { - return 3 + TMath::Min(m1, m2); + return 3 + std::min(m1, m2); } if (d == 2) { - return 6 + TMath::Min(m1, m2); + return 6 + std::min(m1, m2); } return 9; } //_____________________________________________________________________________ - bool TestLambda(float pt, float l1, float l2) + bool testLambda(float pt, float l1, float l2) { // Parameterization for full dispersion // Parameterizatino for full dispersion float l2Mean = 1.53126 + 9.50835e+06 / (1. + 1.08728e+07 * pt + 1.73420e+06 * pt * pt); - float l1Mean = 1.12365 + 0.123770 * TMath::Exp(-pt * 0.246551) + 5.30000e-03 * pt; + float l1Mean = 1.12365 + 0.123770 * std::exp(-pt * 0.246551) + 5.30000e-03 * pt; float l2Sigma = 6.48260e-02 + 7.60261e+10 / (1. + 1.53012e+11 * pt + 5.01265e+05 * pt * pt) + 9.00000e-03 * pt; float l1Sigma = 4.44719e-04 + 6.99839e-01 / (1. + 1.22497e+00 * pt + 6.78604e-07 * pt * pt) + 9.00000e-03 * pt; - float c = -0.35 - 0.550 * TMath::Exp(-0.390730 * pt); + float c = -0.35 - 0.550 * std::exp(-0.390730 * pt); return 0.5 * (l1 - l1Mean) * (l1 - l1Mean) / l1Sigma / l1Sigma + 0.5 * (l2 - l2Mean) * (l2 - l2Mean) / l2Sigma / l2Sigma + @@ -621,8 +704,8 @@ struct phosPi0 { if (res < 0) return 0; - if (res >= nMaxMixBins) - return nMaxMixBins - 1; + if (res >= kMaxMixBins) + return kMaxMixBins - 1; return res; } //---------------------------------------- @@ -639,13 +722,13 @@ struct phosPi0 { return mcParticles->iteratorAt(iparent1).pdgCode(); } auto parent2 = mcParticles->iteratorAt(iparent2); - if (parent2.mothersIds().size() == 0 || parent2.pdgCode() == 21 || abs(parent2.pdgCode()) < 11 || abs(parent2.pdgCode()) > 5000) { // no parents, parent not quark/gluon, strings + if (parent2.mothersIds().size() == 0 || parent2.pdgCode() == 21 || std::abs(parent2.pdgCode()) < 11 || std::abs(parent2.pdgCode()) > 5000) { // no parents, parent not quark/gluon, strings break; } iparent2 = parent2.mothersIds()[0]; } auto parent1 = mcParticles->iteratorAt(iparent1); - if (parent1.mothersIds().size() == 0 || parent1.pdgCode() == 21 || abs(parent1.pdgCode()) < 11 || abs(parent1.pdgCode()) > 5000) { // no parents, parent not quark/gluon, strings + if (parent1.mothersIds().size() == 0 || parent1.pdgCode() == 21 || std::abs(parent1.pdgCode()) < 11 || std::abs(parent1.pdgCode()) > 5000) { // no parents, parent not quark/gluon, strings return 0; } iparent1 = parent1.mothersIds()[0]; From aa07a72d759b98ea43cc8a2724b6ef6a5ccd929f Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 17 Jan 2025 18:34:57 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- Common/TableProducer/caloClusterProducer.cxx | 30 ++++++++++---------- PWGEM/Tasks/phosNonlin.cxx | 2 +- PWGEM/Tasks/phosPi0.cxx | 22 +++++++------- 3 files changed, 27 insertions(+), 27 deletions(-) diff --git a/Common/TableProducer/caloClusterProducer.cxx b/Common/TableProducer/caloClusterProducer.cxx index a6afbd901db..9516458f1f0 100644 --- a/Common/TableProducer/caloClusterProducer.cxx +++ b/Common/TableProducer/caloClusterProducer.cxx @@ -76,8 +76,8 @@ struct caloClusterProducerTask { static constexpr int16_t kCpvX = 7; // grid 6 steps along z and 7 along phi as largest match ellips 20x20 cm static constexpr int16_t kCpvZ = 6; static constexpr int16_t kCpvCells = 4 * kCpvX * kCpvZ; // 4 modules - static constexpr float kCpvMaxX = 73; // max CPV coordinate phi - static constexpr float kCpvMaxZ = 63; // max CPV coordinate z + static constexpr float kCpvMaxX = 73; // max CPV coordinate phi + static constexpr float kCpvMaxZ = 63; // max CPV coordinate z class TrackMatch { @@ -128,7 +128,7 @@ struct caloClusterProducerTask { } std::map bcMap; int bcId = 0; - for (auto & bc : bcs) { + for (auto& bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -136,7 +136,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (const auto & cl : colls) { + for (const auto& cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -351,7 +351,7 @@ struct caloClusterProducerTask { float sigmaX = 1. / std::min(5.2, 1.111 + 0.56 * std::exp(-0.031 * e * e) + 4.8 / std::pow(e + 0.61, 3)); // inverse sigma X float sigmaZ = 1. / std::min(3.3, 1.12 + 0.35 * std::exp(-0.032 * e * e) + 0.75 / std::pow(e + 0.24, 3)); // inverse sigma Z - for (const int & indx : regions) { + for (const int& indx : regions) { if (indx >= 0 && indx < kCpvCells) { for (int ii = cpvPoints->mStart[indx]; ii < cpvPoints->mEnd[indx]; ii++) { auto p = cpvMatchPoints[indx][ii]; @@ -422,7 +422,7 @@ struct caloClusterProducerTask { } std::map bcMap; int bcId = 0; - for (auto const & bc : bcs) { + for (auto const& bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -430,7 +430,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (auto const &cl : colls) { + for (auto const& cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -694,7 +694,7 @@ struct caloClusterProducerTask { mclabels.clear(); mcamplitudes.clear(); gsl::span spDigList = outputTruthCont.getLabels(i); - for (const auto &cellLab : spDigList) { + for (const auto& cellLab : spDigList) { mclabels.push_back(cellLab.getTrackID()); // Track ID in current event? mcamplitudes.push_back(cellLab.getEdep()); } @@ -765,7 +765,7 @@ struct caloClusterProducerTask { std::map bcMap; int bcId = 0; - for (const auto &bc : bcs) { + for (const auto& bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -773,7 +773,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (const auto &cl : colls) { + for (const auto& cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -1174,7 +1174,7 @@ struct caloClusterProducerTask { std::map bcMap; int bcId = 0; - for (const auto & bc : bcs) { + for (const auto& bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -1182,7 +1182,7 @@ struct caloClusterProducerTask { // If several collisions appear in BC, choose one with largers number of contributors std::map colMap; int colId = 0; - for (const auto & cl : colls) { + for (const auto& cl : colls) { auto colbc = colMap.find(cl.bc_as().globalBC()); if (colbc == colMap.end()) { // single collision per BC colMap[cl.bc_as().globalBC()] = colId; @@ -1545,7 +1545,7 @@ struct caloClusterProducerTask { mclabels.clear(); mcamplitudes.clear(); gsl::span spDigList = outputTruthCont.getLabels(i); - for (const auto & cellLab : spDigList) { + for (const auto& cellLab : spDigList) { mclabels.push_back(cellLab.getTrackID()); // Track ID in current event? mcamplitudes.push_back(cellLab.getEdep()); } @@ -1617,7 +1617,7 @@ struct caloClusterProducerTask { const float etaMax = 0.178266; double bz = o2::base::Propagator::Instance()->getNominalBz(); // magnetic field - trackPhi = RecoDecay::constrainAngle(trackPhi,0.,1); //constrain angle to range 0,twoPi + trackPhi = RecoDecay::constrainAngle(trackPhi, 0., 1); // constrain angle to range 0,twoPi if (trackPhi < phiMin || trackPhi > phiMax || std::abs(trackEta) > etaMax) { // do not match even approximately return false; } @@ -1632,7 +1632,7 @@ struct caloClusterProducerTask { } // get PHOS radius - const double shiftY = -1.26; // Depth-optimized + const double shiftY = -1.26; // Depth-optimized double posL[3] = {0., 0., shiftY}; // local position at the center of module double posG[3] = {0}; geomPHOS->getAlignmentMatrix(module)->LocalToMaster(posL, posG); diff --git a/PWGEM/Tasks/phosNonlin.cxx b/PWGEM/Tasks/phosNonlin.cxx index a0d28285326..edaca3db5c3 100644 --- a/PWGEM/Tasks/phosNonlin.cxx +++ b/PWGEM/Tasks/phosNonlin.cxx @@ -180,7 +180,7 @@ struct phosNonlin { mCurEvent.clear(); int i, j, k, l; - for (auto const & clu : clusters) { + for (auto const& clu : clusters) { if (clu.e() < minCluE || clu.ncell() < minCluNcell || clu.time() > maxCluTime || clu.time() < minCluTime || diff --git a/PWGEM/Tasks/phosPi0.cxx b/PWGEM/Tasks/phosPi0.cxx index 77aeb76f081..8d0c1cf39b4 100644 --- a/PWGEM/Tasks/phosPi0.cxx +++ b/PWGEM/Tasks/phosPi0.cxx @@ -374,7 +374,7 @@ struct phosPi0 { if constexpr (isMC) { // check current collision Id for clusters int cluMcBCId = -1; - for (const auto &clu : clusters) { + for (const auto& clu : clusters) { auto mcList = clu.labels(); // const std::vector int nParents = mcList.size(); for (int iParent = 0; iParent < nParents; iParent++) { // Not found nbar parent yiet @@ -462,7 +462,7 @@ struct phosPi0 { } Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel); // Mix with other photons added to stack - for (cosnt auto & ph2 : mCurEvent) { + for (cosnt auto& ph2 : mCurEvent) { double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { @@ -510,9 +510,9 @@ struct phosPi0 { } // Mixed - for (const auto & ph1 : mCurEvent) { - for (const auto & mixEvent : mMixedEvents[mixedEventBin]) { - for (const auto &ph2 : mixEvent) { + for (const auto& ph1 : mCurEvent) { + for (const auto& mixEvent : mMixedEvents[mixedEventBin]) { + for (const auto& ph2 : mixEvent) { double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { @@ -597,14 +597,14 @@ struct phosPi0 { int mcLabel = -1; Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel); // Mix with other photons added to stack - for (const auto & ph2 : mCurEvent) { + for (const auto& ph2 : mCurEvent) { double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { m = std::sqrt(m); } double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) + - std::pow(ph1.py + ph2.py, 2)); + std::pow(ph1.py + ph2.py, 2)); int modComb = moduleCombination(ph1.mod, ph2.mod); hReMod->Fill(m, pt, modComb); hReAll->Fill(m, pt); @@ -625,16 +625,16 @@ struct phosPi0 { } // Mixed - for (const auto & ph1 : mCurEvent) { - for (const auto & mixEvent : mMixedEvents[mixedEventBin]) { - for (const auto &ph2 : mixEvent) { + for (const auto& ph1 : mCurEvent) { + for (const auto& mixEvent : mMixedEvents[mixedEventBin]) { + for (const auto& ph2 : mixEvent) { double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) - std::pow(ph1.py + ph2.py, 2) - std::pow(ph1.pz + ph2.pz, 2); if (m > 0) { m = std::sqrt(m); } double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) + - std::pow(ph1.py + ph2.py, 2)); + std::pow(ph1.py + ph2.py, 2)); int modComb = moduleCombination(ph1.mod, ph2.mod); hMiMod->Fill(m, pt, modComb); hMiAll->Fill(m, pt);