[PWGJE] modification for systematic uncertainty of b-jets and multiplicity analysis (#15830)

Author: hanseopark

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 <CCDB/BasicCCDBManager.h>
-#include <CommonConstants/PhysicsConstants.h>
-#include <DCAFitter/DCAFitterN.h>
 #include <DetectorsBase/MatLayerCylSet.h>
 #include <DetectorsBase/Propagator.h>
-#include <Framework/ASoA.h>
-#include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
-#include <Framework/AnalysisTask.h>
 #include <Framework/Configurable.h>
 #include <Framework/HistogramRegistry.h>
 #include <Framework/HistogramSpec.h>
 #include <Framework/InitContext.h>
 #include <Framework/Logger.h>
 #include <Framework/runDataProcessing.h>
-#include <ReconstructionDataFormats/DCA.h>
 #include <ReconstructionDataFormats/TrackParametrizationWithError.h>
 #include <ReconstructionDataFormats/Vertex.h>
 
@@ -74,12 +74,12 @@ struct SecondaryVertexReconstruction {
   Configurable<bool> propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"};
   Configurable<bool> useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"};
   Configurable<bool> useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"};
-  Configurable<double> maxR{"maxR", 200., "reject PCA's above this radius"};
-  Configurable<double> maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"};
-  Configurable<double> maxRsv{"maxRsv", 999., "max. radius of the reconstruced SV"};
-  Configurable<double> maxZsv{"maxZsv", 999., "max. Z coordinates of the reconstruced SV"};
-  Configurable<double> minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"};
-  Configurable<double> minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"};
+  Configurable<float> maxR{"maxR", 200., "reject PCA's above this radius"};
+  Configurable<float> maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"};
+  Configurable<float> maxRsv{"maxRsv", 999., "max. radius of the reconstruced SV"};
+  Configurable<float> maxZsv{"maxZsv", 999., "max. Z coordinates of the reconstruced SV"};
+  Configurable<float> minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"};
+  Configurable<float> minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"};
   Configurable<float> ptMinTrack{"ptMinTrack", -1., "min. track pT"};
   Configurable<float> etaMinTrack{"etaMinTrack", -99999., "min. pseudorapidity"};
   Configurable<float> 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<o2::ccdb::BasicCCDBManager> 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<std::array<float, 3>, numProngs> arrayMomenta;
+      std::array<std::array<float, 3>, numProngs> arrayMomenta{};
       std::array<o2::dataformats::DCA, numProngs> 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<double, numProngs> massArray;
+      std::array<double, numProngs> 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<SecondaryVertexReconstruction>(cfgc, TaskName{"jet-sv-reconstruction-charged"})}; // o2-linter: disable=name/o2-task
+  return WorkflowSpec{adaptAnalysisTask<SecondaryVertexReconstruction>(cfgc, TaskName{"jet-sv-reconstruction-charged"})}; // o2-linter: disable=name/o2-task (templated struct)
 }

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 <typename T, typename U>
+  void fillRhoAreaSubtractedHistogramSV3ProngData(T const& jet, U const& /*prongs*/, float centrality, float rho)
+  {
+    if (jet.template secondaryVertices_as<U>().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<U>().size(), centrality);
+      for (const auto& prong : jet.template secondaryVertices_as<U>()) {
+        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<U>(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<U>(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 <typename T, typename U>
   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<soa::Join<aod::JetCollisions, aod::BkgChargedRhos>>::iterator const& collision, soa::Join<JetTableData, TagTableData, aod::DataSecondaryVertex3ProngIndices> 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<aod::JetTracks>(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<soa::Join<aod::JetCollisionsMCD, aod::JCollisionOutliers>>::iterator const& collision, soa::Join<JetTableMCD, TagTableMCD, aod::MCDSecondaryVertex3ProngIndices> const& mcdjets, aod::MCDSecondaryVertex3Prongs const& prongs)
   {
     if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {

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<JetTagTracksData>()) {
         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<aod::JetTracks>(mcdjet)) {
         continue;
       }
-      auto eventWeight = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>().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<aod::JetTracks>(mcdjet)) {
         continue;
       }
-      auto eventWeight = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>().weight();
-      fillValidationFlavourDefMCD<soa::Join<JetTableMCP, JetTableMCPMCD>>(mcdjet, tracks, particles, particlesPerColl, eventWeight);
+      fillValidationFlavourDefMCD<soa::Join<JetTableMCP, JetTableMCPMCD>>(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);