diff --git a/Common/TableProducer/caloClusterProducer.cxx b/Common/TableProducer/caloClusterProducer.cxx index b569cf7a58f..cbebf80a53b 100644 --- a/Common/TableProducer/caloClusterProducer.cxx +++ b/Common/TableProducer/caloClusterProducer.cxx @@ -9,6 +9,17 @@ // 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 +#include +#include +#include +#include + #include "Framework/ConfigParamSpec.h" #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" @@ -36,9 +47,9 @@ using namespace o2::framework; using namespace o2; -using mcCells = o2::soa::Join; +using McCells = o2::soa::Join; -struct caloClusterProducerTask { +struct CaloClusterProducer { Produces clucursor; Produces cluambcursor; Produces matchedTracks; @@ -46,13 +57,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 +82,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 +98,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 +134,7 @@ struct caloClusterProducerTask { } std::map bcMap; int bcId = 0; - for (auto bc : bcs) { + for (const auto& bc : bcs) { bcMap[bc.globalBC()] = bcId; bcId++; } @@ -131,7 +142,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 +160,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 +193,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 +227,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 +256,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 +266,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 +304,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 +317,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; } @@ -358,8 +369,8 @@ struct caloClusterProducerTask { } } } - if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + if (cpvdist != 99.) { // was evaluated + cpvdist = std::sqrt(cpvdist); // was squared } int cpvindex = -2; // -2 no CPV in event if (cpvExist) { @@ -400,11 +411,11 @@ struct caloClusterProducerTask { } } - PROCESS_SWITCH(caloClusterProducerTask, processStandalone, "Process PHOS and CPV only", true); + PROCESS_SWITCH(CaloClusterProducer, processStandalone, "Process PHOS and CPV only", true); 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 +428,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 +436,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 +454,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 +488,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 +543,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 +572,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 +582,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 +620,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 +633,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; } @@ -674,8 +685,8 @@ struct caloClusterProducerTask { } } } - if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + if (cpvdist != 99.) { // was evaluated + cpvdist = std::sqrt(cpvdist); // was squared } int cpvindex = -2; // -2 no CPV in event if (cpvExist) { @@ -689,7 +700,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()); } @@ -730,7 +741,7 @@ struct caloClusterProducerTask { } } - PROCESS_SWITCH(caloClusterProducerTask, processStandaloneMC, "Process MC, PHOS and CPV only", true); + PROCESS_SWITCH(CaloClusterProducer, processStandaloneMC, "Process MC, PHOS and CPV only", false); //------------------------------------------------------------ void processFull(o2::aod::BCsWithTimestamps const& bcs, @@ -760,7 +771,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 +779,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 +796,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 +829,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 +863,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 +892,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 +901,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 +914,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 +941,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 +965,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 +976,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 +1023,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 +1036,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 +1094,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; @@ -1094,11 +1105,11 @@ struct caloClusterProducerTask { } } - if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + if (cpvdist != 99.) { // was evaluated + cpvdist = std::sqrt(cpvdist); // was squared } - if (trackdist != 99.) { // was evaluated - trackdist = sqrt(trackdist); // was squared + if (trackdist != 99.) { // was evaluated + trackdist = std::sqrt(trackdist); // was squared } float lambdaShort = 0., lambdaLong = 0.; @@ -1140,12 +1151,12 @@ struct caloClusterProducerTask { } } - PROCESS_SWITCH(caloClusterProducerTask, processFull, "Process with track matching", false); + PROCESS_SWITCH(CaloClusterProducer, processFull, "Process with track matching", false); //------------------------------------------------------------ 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 +1180,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 +1188,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 +1205,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 +1239,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 +1294,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 +1323,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 +1332,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 +1345,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 +1372,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 +1396,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 +1407,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 +1453,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 +1466,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 +1522,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; @@ -1522,11 +1533,11 @@ struct caloClusterProducerTask { } } - if (cpvdist != 99.) { // was evaluated - cpvdist = sqrt(cpvdist); // was squared + if (cpvdist != 99.) { // was evaluated + cpvdist = std::sqrt(cpvdist); // was squared } - if (trackdist != 99.) { // was evaluated - trackdist = sqrt(trackdist); // was squared + if (trackdist != 99.) { // was evaluated + trackdist = std::sqrt(trackdist); // was squared } float lambdaShort = 0., lambdaLong = 0.; @@ -1540,7 +1551,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()); } @@ -1581,17 +1592,17 @@ struct caloClusterProducerTask { } } - PROCESS_SWITCH(caloClusterProducerTask, processFullMC, "Process MC with track matching", false); + PROCESS_SWITCH(CaloClusterProducer, 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 +1623,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 +1638,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 +1667,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 +1676,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 +1693,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 +1709,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 +1721,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(static_cast(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(static_cast(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; @@ -1743,5 +1773,5 @@ struct caloClusterProducerTask { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc)}; } diff --git a/PWGEM/Tasks/CMakeLists.txt b/PWGEM/Tasks/CMakeLists.txt index 4f0a92bd6f2..ef263b02e2e 100644 --- a/PWGEM/Tasks/CMakeLists.txt +++ b/PWGEM/Tasks/CMakeLists.txt @@ -9,47 +9,47 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. -o2physics_add_dpl_workflow(phoscellqa +o2physics_add_dpl_workflow(phos-cell-q-a SOURCES phosCellQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(phoscluqa +o2physics_add_dpl_workflow(phos-clu-q-a SOURCES phosCluQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(phostrigqa +o2physics_add_dpl_workflow(phos-trig-q-a SOURCES phosTrigQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(phospi0 +o2physics_add_dpl_workflow(phos-pi0 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 +o2physics_add_dpl_workflow(phos-calibration SOURCES phosCalibration.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase O2::PHOSReconstruction COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(phosalign +o2physics_add_dpl_workflow(phos-align SOURCES phosAlign.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase O2::PHOSReconstruction COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(phosnbar +o2physics_add_dpl_workflow(phos-nbar SOURCES phosNbar.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(phosnonlin +o2physics_add_dpl_workflow(phos-nonlin 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 +o2physics_add_dpl_workflow(phos-el-id SOURCES phosElId.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::PHOSBase COMPONENT_NAME Analysis) diff --git a/PWGEM/Tasks/phosNonlin.cxx b/PWGEM/Tasks/phosNonlin.cxx index 9115fa432c8..f25043c5f0c 100644 --- a/PWGEM/Tasks/phosNonlin.cxx +++ b/PWGEM/Tasks/phosNonlin.cxx @@ -9,17 +9,26 @@ // 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 #include #include +#include #include -#include + #include #include "Common/DataModel/CaloClusters.h" #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 +43,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"}; +struct PhosNonlin { + Configurable skimmedProcessing{"skimmedProcessing", true, "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", 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,29 +87,34 @@ 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; - std::array hReIJ, hReKL, hReMIJ, hReMKL; + static constexpr int kNp = 10; + std::array hReIJ, hReKL, hReMIJ, hReMKL; TH2* hMi; std::vector pt = {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, @@ -116,21 +126,26 @@ struct phosNonlin { { LOG(info) << "Initializing PHOS nonlin analysis task ..."; - const AxisSpec - ptAxis{pt, "p_{T} (GeV/c)"}, + 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})"}; - for (int i = 0; i < mNp; i++) { - for (int j = 0; j < mNp; j++) { - hReIJ[i * mNp + j] = std::get>(mHistManager1.add(Form("hRe_a%d_b%d", i, j), "inv mass", + for (int i = 0; i < kNp; i++) { + for (int j = 0; j < kNp; j++) { + hReIJ[i * kNp + j] = std::get>(mHistManager1.add(Form("hRe_a%d_b%d", i, j), "inv mass", HistType::kTH2F, {mggAxis, ptAxis})) .get(); - hReKL[i * mNp + j] = std::get>(mHistManager1.add(Form("hRe_c%d_d%d", i, j), "inv mass", + hReKL[i * kNp + j] = std::get>(mHistManager1.add(Form("hRe_c%d_d%d", i, j), "inv mass", HistType::kTH2F, {mggAxis, ptAxis})) .get(); - // hReMIJ[i*mNp+j] = std::get>(mHistManager2.add(Form("hReM_a%d_b%d",i,j), "inv mass", + // hReMIJ[i*kNp+j] = std::get>(mHistManager2.add(Form("hReM_a%d_b%d",i,j), "inv mass", // HistType::kTH2F, {mggAxis, ptAxis})).get(); - // hReMKL[i*mNp+j] = std::get>(mHistManager2.add(Form("hReM_c%d_d%d",i,j), "inv mass", + // hReMKL[i*kNp+j] = std::get>(mHistManager2.add(Form("hReM_c%d_d%d",i,j), "inv mass", // HistType::kTH2F, {mggAxis, ptAxis})).get(); } } @@ -141,57 +156,71 @@ 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(); - k = mNp / 2; - 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)); - hReIJ[i * mNp + j]->Fill(m12, pt1); - hReIJ[i * mNp + j]->Fill(m12, pt2); + k = kNp / 2; + l = kNp / 2; + for (i = 0; i < kNp; i++) { + for (j = 0; j < kNp; j++) { + 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 * kNp + j]->Fill(m12, pt1); + hReIJ[i * kNp + j]->Fill(m12, pt2); // if(ph1.mod==ph2.mod){ - // hReMIJ[i*mNp + j]->Fill(m12,pt1); - // hReMIJ[i*mNp + j]->Fill(m12,pt2); + // hReMIJ[i*kNp + j]->Fill(m12,pt1); + // hReMIJ[i*kNp + j]->Fill(m12,pt2); // } } } } - i = mNp / 2; - 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)); - hReKL[k * mNp + l]->Fill(m12, pt1); - hReKL[k * mNp + l]->Fill(m12, pt2); + i = kNp / 2; + j = kNp / 2; + for (k = 0; k < kNp; k++) { + for (l = 0; l < kNp; 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 * kNp + l]->Fill(m12, pt1); + hReKL[k * kNp + l]->Fill(m12, pt2); // if(ph1.mod==ph2.mod){ - // hReMKL[k*mNp + l]->Fill(m12,pt1); - // hReMKL[k*mNp + l]->Fill(m12,pt2); + // hReMKL[k*kNp + l]->Fill(m12,pt1); + // hReMKL[k*kNp + l]->Fill(m12,pt2); // } } } @@ -202,19 +231,19 @@ struct phosNonlin { } // Mixed - for (auto ph1 : mCurEvent) { - for (auto mixEvent : mMixedEvents[mixedEventBin]) { - for (auto ph2 : mixEvent) { + for (const auto& ph1 : mCurEvent) { + for (const auto& mixEvent : mMixedEvents[mixedEventBin]) { + for (const auto& ph2 : mixEvent) { double m = (ph1 + ph2).M(); double pt1 = ph1.Pt(); double pt2 = ph2.Pt(); - i = mNp / 2; - 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)); + i = kNp / 2; + j = kNp / 2; + k = kNp / 2; + l = kNp / 2; + 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 +254,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 +269,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 - kNp / 2) + mdAj * (j - kNp / 2); + const double b = mB + mdBi * (i - kNp / 2) + mdBj * (j - kNp / 2); + const double c = mC + mdCi * (i - kNp / 2) + mdCj * (j - kNp / 2); + const double d = mD + mdDk * (k - kNp / 2) + mdDl * (l - kNp / 2); + const double e = mE + mdEk * (k - kNp / 2) + mdEl * (l - kNp / 2); + const double f = mF + mdFk * (k - kNp / 2) + mdFl * (l - kNp / 2); + const double g = mG + mdGk * (k - kNp / 2) + mdGl * (l - kNp / 2); + const double s = mS + mdSi * (i - kNp / 2) + mdSj * (j - kNp / 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 - kNp / 2) + mdAj * (j - kNp / 2); + const double b = mB + mdBi * (i - kNp / 2) + mdBj * (j - kNp / 2); + const double c = mC + mdCi * (i - kNp / 2) + mdCj * (j - kNp / 2); + const double d = mD + mdDk * (k - kNp / 2) + mdDl * (l - kNp / 2); + const double e = mE + mdEk * (k - kNp / 2) + mdEl * (l - kNp / 2); + return a + b / std::sqrt(en) + c / en + d / (en * std::sqrt(en)) + e / (en * en); } return 0.; } @@ -275,5 +304,5 @@ struct phosNonlin { o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc) { return o2::framework::WorkflowSpec{ - o2::framework::adaptAnalysisTask(cfgc)}; + o2::framework::adaptAnalysisTask(cfgc)}; } diff --git a/PWGEM/Tasks/phosPi0.cxx b/PWGEM/Tasks/phosPi0.cxx index fe248dbeaa6..529a1f29284 100644 --- a/PWGEM/Tasks/phosPi0.cxx +++ b/PWGEM/Tasks/phosPi0.cxx @@ -9,11 +9,19 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file phosPi0.cxx +/// \brief PHOS pi0/eta analysis +/// \author Dmitri Peresunko +/// + +#include #include #include #include #include +#include #include + #include "Common/DataModel/CaloClusters.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" @@ -27,56 +35,57 @@ #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"}; +struct PhosPi0 { + 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"}; + 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; } + bool isCPVOK() const { return (mPID >> 2) & 1; } + bool isDispOK() const { return (mPID >> 1) & 1; } public: double px = 0.; // px @@ -88,11 +97,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 +120,14 @@ 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)"}, + 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 +138,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 +162,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 +195,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,74 +228,147 @@ 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::math::TwoPI, "#phi (rad)"}}); + mHistManager.add("hMCPi0SecVtx", "pi0 secondary", HistType::kTH2F, {{100, 0., 500., "R (cm)"}, {100, -o2::constants::math::PI, o2::constants::math::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); + 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); } - PROCESS_SWITCH(phosPi0, processMC, "processMC", false); + PROCESS_SWITCH(PhosPi0, processMC, "processMC", false); template void scanAll(TCollision& col, 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) && std::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 +376,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 +394,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 +424,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 +462,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 (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); bool isPi0 = false; @@ -427,17 +512,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 +541,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 +565,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 +597,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 +627,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,41 +656,41 @@ 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(); } } } - PROCESS_SWITCH(phosPi0, processBC, "processBC", false); + 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 +706,8 @@ struct phosPi0 { if (res < 0) return 0; - if (res >= nMaxMixBins) - return nMaxMixBins - 1; + if (res >= kMaxMixBins) + return kMaxMixBins - 1; return res; } //---------------------------------------- @@ -639,13 +724,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]; @@ -657,5 +742,5 @@ struct phosPi0 { o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc) { return o2::framework::WorkflowSpec{ - o2::framework::adaptAnalysisTask(cfgc)}; + o2::framework::adaptAnalysisTask(cfgc)}; }