diff --git a/PWGJE/TableProducer/secondaryVertexReconstruction.cxx b/PWGJE/TableProducer/secondaryVertexReconstruction.cxx index e2f2246b21b..1afe05455f2 100644 --- a/PWGJE/TableProducer/secondaryVertexReconstruction.cxx +++ b/PWGJE/TableProducer/secondaryVertexReconstruction.cxx @@ -22,22 +22,22 @@ #include "Common/Core/trackUtilities.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "Framework/ASoA.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "ReconstructionDataFormats/DCA.h" #include -#include -#include #include #include -#include -#include #include -#include #include #include #include #include #include #include -#include #include #include @@ -74,12 +74,12 @@ struct SecondaryVertexReconstruction { Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; - Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; - Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; - Configurable maxRsv{"maxRsv", 999., "max. radius of the reconstruced SV"}; - Configurable maxZsv{"maxZsv", 999., "max. Z coordinates of the reconstruced SV"}; - Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; - Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; + Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable maxRsv{"maxRsv", 999., "max. radius of the reconstruced SV"}; + Configurable maxZsv{"maxZsv", 999., "max. Z coordinates of the reconstruced SV"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; Configurable ptMinTrack{"ptMinTrack", -1., "min. track pT"}; Configurable etaMinTrack{"etaMinTrack", -99999., "min. pseudorapidity"}; Configurable etaMaxTrack{"etaMaxTrack", 4., "max. pseudorapidity"}; @@ -94,13 +94,15 @@ struct SecondaryVertexReconstruction { o2::vertexing::DCAFitterN<2> df2; // 2-prong vertex fitter o2::vertexing::DCAFitterN<3> df3; // 3-prong vertex fitter Service ccdb; - o2::base::MatLayerCylSet* lut; + o2::base::MatLayerCylSet* lut = nullptr; HistogramRegistry registry{"registry"}; int runNumber{0}; float toMicrometers = 10000.; // from cm to µm double bz{0.}; + static constexpr int TwoProngCount = 2; + static constexpr int ThreeProngCount = 3; void init(InitContext const&) { @@ -216,7 +218,7 @@ struct SecondaryVertexReconstruction { auto covMatrixPV = primaryVertex.getCov(); // Get track momenta and impact parameters - std::array, numProngs> arrayMomenta; + std::array, numProngs> arrayMomenta{}; std::array impactParameters; for (unsigned int inum = 0; inum < numProngs; ++inum) { trackParVars[inum].getPxPyPzGlo(arrayMomenta[inum]); @@ -236,12 +238,12 @@ struct SecondaryVertexReconstruction { auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); // calculate invariant mass - std::array massArray; + std::array massArray{}; std::fill(massArray.begin(), massArray.end(), o2::constants::physics::MassPiPlus); - double massSV = RecoDecay::m(std::move(arrayMomenta), massArray); + double massSV = RecoDecay::m(arrayMomenta, massArray); // fill candidate table rows - if ((doprocessData3Prongs || doprocessData3ProngsExternalMagneticField) && numProngs == 3) { + if ((doprocessData3Prongs || doprocessData3ProngsExternalMagneticField) && numProngs == ThreeProngCount) { sv3prongTableData(analysisJet.globalIndex(), primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], @@ -250,7 +252,7 @@ struct SecondaryVertexReconstruction { arrayMomenta[0][2] + arrayMomenta[1][2] + arrayMomenta[2][2], energySV, massSV, chi2PCA, dispersion, errorDecayLength, errorDecayLengthXY); svIndices.push_back(sv3prongTableData.lastIndex()); - } else if ((doprocessData2Prongs || doprocessData2ProngsExternalMagneticField) && numProngs == 2) { + } else if ((doprocessData2Prongs || doprocessData2ProngsExternalMagneticField) && numProngs == TwoProngCount) { sv2prongTableData(analysisJet.globalIndex(), primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], @@ -259,7 +261,7 @@ struct SecondaryVertexReconstruction { arrayMomenta[0][2] + arrayMomenta[1][2], energySV, massSV, chi2PCA, dispersion, errorDecayLength, errorDecayLengthXY); svIndices.push_back(sv2prongTableData.lastIndex()); - } else if ((doprocessMCD3Prongs || doprocessMCD3ProngsExternalMagneticField) && numProngs == 3) { + } else if ((doprocessMCD3Prongs || doprocessMCD3ProngsExternalMagneticField) && numProngs == ThreeProngCount) { sv3prongTableMCD(analysisJet.globalIndex(), primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], @@ -268,7 +270,7 @@ struct SecondaryVertexReconstruction { arrayMomenta[0][2] + arrayMomenta[1][2] + arrayMomenta[2][2], energySV, massSV, chi2PCA, dispersion, errorDecayLength, errorDecayLengthXY); svIndices.push_back(sv3prongTableMCD.lastIndex()); - } else if ((doprocessMCD2Prongs || doprocessMCD2ProngsExternalMagneticField) && numProngs == 2) { + } else if ((doprocessMCD2Prongs || doprocessMCD2ProngsExternalMagneticField) && numProngs == TwoProngCount) { sv2prongTableMCD(analysisJet.globalIndex(), primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], @@ -399,5 +401,5 @@ struct SecondaryVertexReconstruction { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-sv-reconstruction-charged"})}; // o2-linter: disable=name/o2-task + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"jet-sv-reconstruction-charged"})}; // o2-linter: disable=name/o2-task (templated struct) } diff --git a/PWGJE/Tasks/bjetCentMult.cxx b/PWGJE/Tasks/bjetCentMult.cxx index 576afeb4ac0..8bd8111b200 100644 --- a/PWGJE/Tasks/bjetCentMult.cxx +++ b/PWGJE/Tasks/bjetCentMult.cxx @@ -162,6 +162,25 @@ struct BjetCentMultTask { registry.add("hn_taggedjet_3prong_Sxyz_N1_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxyz}, {axisMass}, {axisCentrality}}}); } } + if (doprocessRhoAreaSubSV3ProngData) { + registry.add("h_event_centrality", "", {HistType::kTH1F, {{axisCentrality}}}); + registry.add("h2_jet_pt_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisJetPt}, {axisCentrality}}}); + registry.add("h2_jet_eta_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisEta}, {axisCentrality}}}); + registry.add("h2_jet_phi_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisPhi}, {axisCentrality}}}); + if (fillGeneralSVQA) { + registry.add("h2_3prong_nprongs_rhoareasubtracted_centrality", "", {HistType::kTH2F, {{axisNprongs}, {axisCentrality}}}); + registry.add("hn_jet_3prong_Sxy_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisLxy}, {axisSigmaLxy}, {axisSxy}, {axisCentrality}}}); + if (fillSVxyz) { + registry.add("hn_jet_3prong_Sxyz_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisLxyz}, {axisSigmaLxyz}, {axisSxyz}, {axisCentrality}}}); + } + } + registry.add("hn_jet_3prong_Sxy_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxy}, {axisMass}, {axisCentrality}}}); + registry.add("hn_taggedjet_3prong_Sxy_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxy}, {axisMass}, {axisCentrality}}}); + if (fillSVxyz) { + registry.add("hn_jet_3prong_Sxyz_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxyz}, {axisMass}, {axisCentrality}}}); + registry.add("hn_taggedjet_3prong_Sxyz_N1_rhoareasubtracted_centrality", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisSxyz}, {axisMass}, {axisCentrality}}}); + } + } if (doprocessSV3ProngMCD || doprocessSV3ProngMCPMCDMatched) { registry.add("h_event_centrality", "", {HistType::kTH1F, {{axisCentrality}}}); registry.add("h3_jet_pt_centrality_flavour", "", {HistType::kTH3F, {{axisJetPt}, {axisCentrality}, {axisJetFlavour}}}); @@ -312,6 +331,47 @@ struct BjetCentMultTask { registry.fill(HIST("h2_jet_phi_part_flavour"), mcpjet.phi(), jetflavour, eventWeight); } + template + void fillRhoAreaSubtractedHistogramSV3ProngData(T const& jet, U const& /*prongs*/, float centrality, float rho) + { + if (jet.template secondaryVertices_as().size() < 1) + return; + registry.fill(HIST("h2_jet_pt_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), centrality); + registry.fill(HIST("h2_jet_eta_rhoareasubtracted_centrality"), jet.eta(), centrality); + registry.fill(HIST("h2_jet_phi_rhoareasubtracted_centrality"), jet.phi(), centrality); + if (fillGeneralSVQA) { + registry.fill(HIST("h2_3prong_nprongs_rhoareasubtracted_centrality"), jet.template secondaryVertices_as().size(), centrality); + for (const auto& prong : jet.template secondaryVertices_as()) { + registry.fill(HIST("hn_jet_3prong_Sxy_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), prong.decayLengthXY(), prong.errorDecayLengthXY(), prong.decayLengthXY() / prong.errorDecayLengthXY(), centrality); + if (fillSVxyz) { + registry.fill(HIST("hn_jet_3prong_Sxyz_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), prong.decayLength(), prong.errorDecayLength(), prong.decayLength() / prong.errorDecayLength(), centrality); + } + } + } + bool checkSv = false; + auto bjetCand = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongCuts->at(0), prongCuts->at(1), prongCuts->at(2), prongCuts->at(4), prongCuts->at(5), false, &checkSv); + if (checkSv && jettaggingutilities::svAcceptance(bjetCand, svDispersionMax)) { + auto maxSxy = bjetCand.decayLengthXY() / bjetCand.errorDecayLengthXY(); + auto massSV = bjetCand.m(); + registry.fill(HIST("hn_jet_3prong_Sxy_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxy, massSV, centrality); + if (jet.isTagged(BJetTaggingMethod::SV)) { + registry.fill(HIST("hn_taggedjet_3prong_Sxy_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxy, massSV, centrality); + } + } + if (fillSVxyz) { + checkSv = false; + auto bjetCandXYZ = jettaggingutilities::jetFromProngMaxDecayLength(jet, prongCuts->at(0), prongCuts->at(1), prongCuts->at(3), prongCuts->at(4), prongCuts->at(5), true, &checkSv); + if (checkSv && jettaggingutilities::svAcceptance(bjetCandXYZ, svDispersionMax)) { + auto maxSxyz = bjetCandXYZ.decayLength() / bjetCandXYZ.errorDecayLength(); + auto massSV = bjetCandXYZ.m(); + registry.fill(HIST("hn_jet_3prong_Sxyz_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxyz, massSV, centrality); + if (jet.isTagged(BJetTaggingMethod::SV3D)) { + registry.fill(HIST("hn_taggedjet_3prong_Sxyz_N1_rhoareasubtracted_centrality"), jet.pt() - (rho * jet.area()), maxSxyz, massSV, centrality); + } + } + } + } + template void fillHistogramSV3ProngData(T const& jet, U const& /*prongs*/, float centrality) { @@ -513,6 +573,26 @@ struct BjetCentMultTask { } PROCESS_SWITCH(BjetCentMultTask, processSV3ProngData, "Fill 3prong imformation for data jets", false); + void processRhoAreaSubSV3ProngData(soa::Filtered>::iterator const& collision, soa::Join const& jets, aod::DataSecondaryVertex3Prongs const& prongs) + { + if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) { + return; + } + float centrality = collision.centFT0M(); + float rho = collision.rho(); + registry.fill(HIST("h_event_centrality"), centrality); + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaCuts->at(0), jetEtaCuts->at(1), trackCuts->at(2), trackCuts->at(3))) { + continue; + } + if (!isAcceptedJet(jet)) { + continue; + } + fillRhoAreaSubtractedHistogramSV3ProngData(jet, prongs, centrality, rho); + } + } + PROCESS_SWITCH(BjetCentMultTask, processRhoAreaSubSV3ProngData, "Fill 3prong imformation for data jets with background subtraction", false); + void processSV3ProngMCD(soa::Filtered>::iterator const& collision, soa::Join const& mcdjets, aod::MCDSecondaryVertex3Prongs const& prongs) { if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) { diff --git a/PWGJE/Tasks/jetTaggerHFQA.cxx b/PWGJE/Tasks/jetTaggerHFQA.cxx index 1dd147c0ff7..dac7777daf7 100644 --- a/PWGJE/Tasks/jetTaggerHFQA.cxx +++ b/PWGJE/Tasks/jetTaggerHFQA.cxx @@ -165,7 +165,7 @@ struct JetTaggerHFQA { } } if (doprocessTracksInJetsData) { - registry.add("h2_track_pt_impact_parameter_xy", "", {HistType::kTH2F, {{axisTrackPt}, {axisImpactParameterXY}}}); + registry.add("h3_jet_pt_track_pt_impact_parameter_xy", "", {HistType::kTH3F, {{axisJetPt}, {axisTrackPt}, {axisImpactParameterXY}}}); } if (doprocessSecondaryContaminationMCD) { registry.add("hn_jet_pt_track_pt_impact_parameter_xy_physical_primary_flavour", "", {HistType::kTHnSparseF, {{axisJetPt}, {axisTrackPt}, {axisImpactParameterXY}, {axisJetFlavour}}}); @@ -1158,7 +1158,7 @@ struct JetTaggerHFQA { } for (auto const& track : jet.template tracks_as()) { float varImpXY = track.dcaXY() * jettaggingutilities::cmTomum; - registry.fill(HIST("h2_track_pt_impact_parameter_xy"), track.pt(), varImpXY); + registry.fill(HIST("h3_jet_pt_track_pt_impact_parameter_xy"), jet.pt(), track.pt(), varImpXY); } } } @@ -1179,7 +1179,7 @@ struct JetTaggerHFQA { if (!isAcceptedJet(mcdjet)) { continue; } - auto eventWeight = collision.mcCollision_as>().weight(); + auto eventWeight = collision.weight(); float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent)); if (mcdjet.pt() > pTHatMaxMCD * pTHat) { continue; @@ -1224,8 +1224,7 @@ struct JetTaggerHFQA { if (!isAcceptedJet(mcdjet)) { continue; } - auto eventWeight = collision.mcCollision_as>().weight(); - fillValidationFlavourDefMCD>(mcdjet, tracks, particles, particlesPerColl, eventWeight); + fillValidationFlavourDefMCD>(mcdjet, tracks, particles, particlesPerColl, collision.weight()); } } PROCESS_SWITCH(JetTaggerHFQA, processValFlavourDefMCD, "to check the validation of jet-flavour definition when compared to distance for mcd jets", false);