PWGLF: add QA histo's + autocorr check (AliceO2Group#4354)

* match default settings to those of the builders * add QA histo's + autocorr check - many QA histograms added to keep track of selection criteria - flipped the definition in TPC PID selection flags: QOL improvement for THnSparse projections - autocorrelation check added for the Omega's, where the bachelor kaon could be "reused" as a proton in the V0 of the other anti-cascade * Please consider the following formatting changes (#8) --------- Co-authored-by: ALICE Builder <[email protected]>
creetz16 · Jan 18, 2024 · 2752084 · 2752084
1 parent 41113f5
commit 2752084
Showing 1 changed file with 125 additions and 13 deletions.
diff --git a/PWGLF/Tasks/cascadecorrelations.cxx b/PWGLF/Tasks/cascadecorrelations.cxx
@@ -54,7 +54,7 @@ using FullTracksExtWithPID = soa::Join<aod::Tracks, aod::TracksExtra, aod::Track
 using FullTracksExtIUWithPID = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TracksDCA, aod::pidTPCPi, aod::pidTPCKa, aod::pidTPCPr>;
 
 // Add a column to the cascdataext table: IsSelected.
-// 0 = not selected, 1 = Xi, 2 = Omega, 3 = both
+// 0 = not selected, 1 = Xi, 2 = both, 3 = Omega
 namespace o2::aod
 {
 namespace cascadeflags
@@ -75,44 +75,119 @@ struct cascadeSelector {
   Configurable<float> tpcNsigmaPion{"tpcNsigmaPion", 3, "TPC NSigma pion <- lambda"};
   Configurable<int> minTPCCrossedRows{"minTPCCrossedRows", 80, "min N TPC crossed rows"}; // TODO: finetune! 80 > 159/2, so no split tracks?
   Configurable<int> minITSClusters{"minITSClusters", 4, "minimum number of ITS clusters"};
-  // Configurable<bool>  doTPConly{"doTPConly", false, "use TPC-only tracks"}; // TODO: maybe do this for high pT only? as cascade decays after IB
 
   // Selection criteria - compatible with core wagon autodetect - copied from cascadeanalysis.cxx
   //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
-  Configurable<double> v0setting_cospa{"v0setting_cospa", 0.95, "v0setting_cospa"};
+  Configurable<double> v0setting_cospa{"v0setting_cospa", 0.995, "v0setting_cospa"};
   Configurable<float> v0setting_dcav0dau{"v0setting_dcav0dau", 1.0, "v0setting_dcav0dau"};
   Configurable<float> v0setting_dcapostopv{"v0setting_dcapostopv", 0.1, "v0setting_dcapostopv"};
   Configurable<float> v0setting_dcanegtopv{"v0setting_dcanegtopv", 0.1, "v0setting_dcanegtopv"};
   Configurable<float> v0setting_radius{"v0setting_radius", 0.9, "v0setting_radius"};
   Configurable<double> cascadesetting_cospa{"cascadesetting_cospa", 0.95, "cascadesetting_cospa"};
   Configurable<float> cascadesetting_dcacascdau{"cascadesetting_dcacascdau", 1.0, "cascadesetting_dcacascdau"};
-  Configurable<float> cascadesetting_dcabachtopv{"cascadesetting_dcabachtopv", 0.1, "cascadesetting_dcabachtopv"};
-  Configurable<float> cascadesetting_cascradius{"cascadesetting_cascradius", 0.5, "cascadesetting_cascradius"};
+  Configurable<float> cascadesetting_dcabachtopv{"cascadesetting_dcabachtopv", 0.05, "cascadesetting_dcabachtopv"};
+  Configurable<float> cascadesetting_cascradius{"cascadesetting_cascradius", 0.9, "cascadesetting_cascradius"};
   Configurable<float> cascadesetting_v0masswindow{"cascadesetting_v0masswindow", 0.01, "cascadesetting_v0masswindow"};
   Configurable<float> cascadesetting_mindcav0topv{"cascadesetting_mindcav0topv", 0.01, "cascadesetting_mindcav0topv"};
   //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
 
+  // TODO: variables as function of Omega mass, only do Xi for now
+  AxisSpec vertexAxis = {200, -10.0f, 10.0f, "cm"};
+  AxisSpec dcaAxis = {100, 0.0f, 10.0f, "cm"};
+  AxisSpec invMassAxis = {1000, 1.0f, 2.0f, "Inv. Mass (GeV/c^{2})"};
+  AxisSpec ptAxis = {100, 0, 15, "#it{p}_{T}"};
+  HistogramRegistry registry{
+    "registry",
+    {
+      // basic selection variables
+      {"hV0Radius", "hV0Radius", {HistType::kTH3F, {{100, 0.0f, 100.0f, "cm"}, invMassAxis, ptAxis}}},
+      {"hCascRadius", "hCascRadius", {HistType::kTH3F, {{100, 0.0f, 100.0f, "cm"}, invMassAxis, ptAxis}}},
+      {"hV0CosPA", "hV0CosPA", {HistType::kTH3F, {{100, 0.95f, 1.0f}, invMassAxis, ptAxis}}},
+      {"hCascCosPA", "hCascCosPA", {HistType::kTH3F, {{100, 0.95f, 1.0f}, invMassAxis, ptAxis}}},
+      {"hDCAPosToPV", "hDCAPosToPV", {HistType::kTH3F, {vertexAxis, invMassAxis, ptAxis}}},
+      {"hDCANegToPV", "hDCANegToPV", {HistType::kTH3F, {vertexAxis, invMassAxis, ptAxis}}},
+      {"hDCABachToPV", "hDCABachToPV", {HistType::kTH3F, {vertexAxis, invMassAxis, ptAxis}}},
+      {"hDCAV0ToPV", "hDCAV0ToPV", {HistType::kTH3F, {vertexAxis, invMassAxis, ptAxis}}},
+      {"hDCAV0Dau", "hDCAV0Dau", {HistType::kTH3F, {dcaAxis, invMassAxis, ptAxis}}},
+      {"hDCACascDau", "hDCACascDau", {HistType::kTH3F, {dcaAxis, invMassAxis, ptAxis}}},
+      {"hLambdaMass", "hLambdaMass", {HistType::kTH3F, {{100, 1.0f, 1.2f, "Inv. Mass (GeV/c^{2})"}, invMassAxis, ptAxis}}},
+
+      // invariant mass per cut, start with Xi
+      {"hMassXi0", "Xi inv mass before selections", {HistType::kTH2F, {invMassAxis, ptAxis}}},
+      {"hMassXi1", "Xi inv mass after TPCnCrossedRows cut", {HistType::kTH2F, {invMassAxis, ptAxis}}},
+      {"hMassXi2", "Xi inv mass after ITSnClusters cut", {HistType::kTH2F, {invMassAxis, ptAxis}}},
+      {"hMassXi3", "Xi inv mass after topo cuts", {HistType::kTH2F, {invMassAxis, ptAxis}}},
+      {"hMassXi4", "Xi inv mass after V0 daughters PID cut", {HistType::kTH2F, {invMassAxis, ptAxis}}},
+      {"hMassXi5", "Xi inv mass after bachelor PID cut", {HistType::kTH2F, {invMassAxis, ptAxis}}},
+
+      // ITS & TPC clusters, with Xi inv mass
+      {"hTPCnCrossedRowsPos", "hTPCnCrossedRowsPos", {HistType::kTH3F, {{160, -0.5, 159.5, "TPC crossed rows"}, invMassAxis, ptAxis}}},
+      {"hTPCnCrossedRowsNeg", "hTPCnCrossedRowsNeg", {HistType::kTH3F, {{160, -0.5, 159.5, "TPC crossed rows"}, invMassAxis, ptAxis}}},
+      {"hTPCnCrossedRowsBach", "hTPCnCrossedRowsBach", {HistType::kTH3F, {{160, -0.5, 159.5, "TPC crossed rows"}, invMassAxis, ptAxis}}},
+      {"hITSnClustersPos", "hITSnClustersPos", {HistType::kTH3F, {{8, -0.5, 7.5, "number of ITS clusters"}, invMassAxis, ptAxis}}},
+      {"hITSnClustersNeg", "hITSnClustersNeg", {HistType::kTH3F, {{8, -0.5, 7.5, "number of ITS clusters"}, invMassAxis, ptAxis}}},
+      {"hITSnClustersBach", "hITSnClustersBach", {HistType::kTH3F, {{8, -0.5, 7.5, "number of ITS clusters"}, invMassAxis, ptAxis}}},
+    },
+  };
+
+  // Keep track of which selections the candidates pass
+  void init(InitContext const&)
+  {
+    auto h = registry.add<TH1>("hSelectionStatus", "hSelectionStatus", HistType::kTH1I, {{10, 0, 10, "status"}});
+    h->GetXaxis()->SetBinLabel(1, "All");
+    h->GetXaxis()->SetBinLabel(2, "nTPC OK");
+    h->GetXaxis()->SetBinLabel(3, "nITS OK");
+    h->GetXaxis()->SetBinLabel(4, "Topo OK");
+    h->GetXaxis()->SetBinLabel(5, "V0 PID OK");
+    h->GetXaxis()->SetBinLabel(6, "Bach PID OK");
+  }
   void process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision, aod::CascDataExt const& Cascades, FullTracksExtIUWithPID const&)
   {
     for (auto& casc : Cascades) {
-      // TODO: make QA histo with info on where cascades fail selections
 
-      // Let's try to do some PID & track quality cuts
       // these are the tracks:
       auto bachTrack = casc.bachelor_as<FullTracksExtIUWithPID>();
       auto posTrack = casc.posTrack_as<FullTracksExtIUWithPID>();
       auto negTrack = casc.negTrack_as<FullTracksExtIUWithPID>();
 
+      // topo variables before cuts:
+      registry.fill(HIST("hV0Radius"), casc.v0radius(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hCascRadius"), casc.cascradius(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hV0CosPA"), casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()), casc.mXi(), casc.pt());
+      registry.fill(HIST("hCascCosPA"), casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()), casc.mXi(), casc.pt());
+      registry.fill(HIST("hDCAPosToPV"), casc.dcapostopv(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hDCANegToPV"), casc.dcanegtopv(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hDCABachToPV"), casc.dcabachtopv(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hDCAV0ToPV"), casc.dcav0topv(collision.posX(), collision.posY(), collision.posZ()), casc.mXi(), casc.pt());
+      registry.fill(HIST("hDCAV0Dau"), casc.dcaV0daughters(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hDCACascDau"), casc.dcacascdaughters(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hLambdaMass"), casc.mLambda(), casc.mXi(), casc.pt());
+
+      registry.fill(HIST("hITSnClustersPos"), posTrack.itsNCls(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hITSnClustersNeg"), negTrack.itsNCls(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hITSnClustersBach"), bachTrack.itsNCls(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hTPCnCrossedRowsPos"), posTrack.tpcNClsCrossedRows(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hTPCnCrossedRowsNeg"), negTrack.tpcNClsCrossedRows(), casc.mXi(), casc.pt());
+      registry.fill(HIST("hTPCnCrossedRowsBach"), bachTrack.tpcNClsCrossedRows(), casc.mXi(), casc.pt());
+
+      registry.fill(HIST("hSelectionStatus"), 0); // all the cascade before selections
+      registry.fill(HIST("hMassXi0"), casc.mXi(), casc.pt());
+
       // TPC N crossed rows
       if (posTrack.tpcNClsCrossedRows() < minTPCCrossedRows || negTrack.tpcNClsCrossedRows() < minTPCCrossedRows || bachTrack.tpcNClsCrossedRows() < minTPCCrossedRows) {
         cascflags(0);
         continue;
       }
+      registry.fill(HIST("hSelectionStatus"), 1); // passes nTPC crossed rows
+      registry.fill(HIST("hMassXi1"), casc.mXi(), casc.pt());
+
       // ITS N clusters
       if (posTrack.itsNCls() < minITSClusters || negTrack.itsNCls() < minITSClusters || bachTrack.itsNCls() < minITSClusters) {
         cascflags(0);
         continue;
       }
+      registry.fill(HIST("hSelectionStatus"), 2); // passes nITS clusters
+      registry.fill(HIST("hMassXi2"), casc.mXi(), casc.pt());
 
       //// TOPO CUTS //// TODO: improve!
       double pvx = collision.posX();
@@ -128,10 +203,12 @@ struct cascadeSelector {
         cascflags(0);
         continue;
       }
+      registry.fill(HIST("hSelectionStatus"), 3); // passes topo
+      registry.fill(HIST("hMassXi3"), casc.mXi(), casc.pt());
 
       // TODO: TOF (for pT > 2 GeV per track?)
 
-      //// TPC ////
+      //// TPC PID ////
       // Lambda check
       if (casc.sign() < 0) {
         // Proton check:
@@ -156,17 +233,25 @@ struct cascadeSelector {
           continue;
         }
       }
+      registry.fill(HIST("hSelectionStatus"), 4); // fails at V0 daughters PID
+      registry.fill(HIST("hMassXi4"), casc.mXi(), casc.pt());
+
       // Bachelor check
       if (TMath::Abs(bachTrack.tpcNSigmaPi()) < tpcNsigmaBachelor) {
         if (TMath::Abs(bachTrack.tpcNSigmaKa()) < tpcNsigmaBachelor) {
           // consistent with both!
-          cascflags(3);
+          cascflags(2);
+          registry.fill(HIST("hSelectionStatus"), 5); // passes bach PID
+          registry.fill(HIST("hMassXi5"), casc.mXi(), casc.pt());
           continue;
         }
         cascflags(1);
+        registry.fill(HIST("hSelectionStatus"), 5); // passes bach PID
+        registry.fill(HIST("hMassXi5"), casc.mXi(), casc.pt());
         continue;
       } else if (TMath::Abs(bachTrack.tpcNSigmaKa()) < tpcNsigmaBachelor) {
-        cascflags(2);
+        cascflags(3);
+        registry.fill(HIST("hSelectionStatus"), 5); // passes bach PID
         continue;
       }
       // if we reach here, the bachelor was neither pion nor kaon
@@ -207,7 +292,8 @@ struct cascadeCorrelations {
       {"hLambdaMass", "hLambdaMass", {HistType::kTH1F, {{1000, 0.0f, 10.0f, "Inv. Mass (GeV/c^{2})"}}}},
 
       {"hSelectionFlag", "hSelectionFlag", {HistType::kTH1I, {selectionFlagAxis}}},
-      {"hAutoCorrelation", "hAutoCorrelation", {HistType::kTH1I, {{4, -0.5f, 3.5f, "Types of autocorrelation"}}}},
+      {"hAutoCorrelation", "hAutoCorrelation", {HistType::kTH1I, {{4, -0.5f, 3.5f, "Types of SS autocorrelation"}}}},
+      {"hAutoCorrelationOS", "hAutoCorrelationOS", {HistType::kTH1I, {{2, -1.f, 1.f, "Charge of OS autocorrelated track"}}}},
       {"hPhi", "hPhi", {HistType::kTH1F, {{100, 0, 2 * PI, "#varphi"}}}},
       {"hEta", "hEta", {HistType::kTH1F, {{100, -2, 2, "#eta"}}}},
 
@@ -237,7 +323,7 @@ struct cascadeCorrelations {
 
     // Some QA on the cascades
     for (auto& casc : Cascades) {
-      if (casc.isSelected() != 2) { // not exclusively an Omega --> consistent with Xi or both
+      if (casc.isSelected() <= 2) { // not exclusively an Omega --> consistent with Xi or both
         if (casc.sign() < 0) {
           registry.fill(HIST("hMassXiMinus"), casc.mXi(), casc.pt());
         } else {
@@ -296,13 +382,39 @@ struct cascadeCorrelations {
 
       // Fill the correct histograms based on same-sign or opposite-sign
       if (trigger.sign() * assoc.sign() < 0) { // opposite-sign
+        // check for autocorrelations between mis-identified kaons (omega bach) and protons (lambda daughter) TODO: improve logic?
+        if (trigger.isSelected() >= 2) {
+          if (trigger.sign() > 0 && trigger.bachelorId() == posIdAssoc) {
+            // K+ from trigger Omega is the same as proton from assoc lambda
+            registry.fill(HIST("hAutoCorrelationOS"), 1);
+            continue;
+          }
+          if (trigger.sign() < 0 && trigger.bachelorId() == negIdAssoc) {
+            // K- from trigger Omega is the same as antiproton from assoc antilambda
+            registry.fill(HIST("hAutoCorrelationOS"), -1);
+            continue;
+          }
+        }
+        if (assoc.isSelected() >= 2) {
+          if (assoc.sign() > 0 && assoc.bachelorId() == posIdTrigg) {
+            // K+ from assoc Omega is the same as proton from trigger lambda
+            registry.fill(HIST("hAutoCorrelationOS"), 1);
+            continue;
+          }
+          if (assoc.sign() < 0 && assoc.bachelorId() == negIdTrigg) {
+            // K- from assoc Omega is the same as antiproton from trigger antilambda
+            registry.fill(HIST("hAutoCorrelationOS"), -1);
+            continue;
+          }
+        }
+
         registry.fill(HIST("hDeltaPhiOS"), dphi);
         registry.fill(HIST("hXiXiOS"), dphi, deta, trigger.pt(), assoc.pt(), invMassXiTrigg, invMassXiAssoc, trigger.isSelected(), assoc.isSelected(), collision.posZ(), collision.multFT0M());
         registry.fill(HIST("hXiOmOS"), dphi, deta, trigger.pt(), assoc.pt(), invMassXiTrigg, invMassOmAssoc, trigger.isSelected(), assoc.isSelected(), collision.posZ(), collision.multFT0M());
         registry.fill(HIST("hOmXiOS"), dphi, deta, trigger.pt(), assoc.pt(), invMassOmTrigg, invMassXiAssoc, trigger.isSelected(), assoc.isSelected(), collision.posZ(), collision.multFT0M());
         registry.fill(HIST("hOmOmOS"), dphi, deta, trigger.pt(), assoc.pt(), invMassOmTrigg, invMassOmAssoc, trigger.isSelected(), assoc.isSelected(), collision.posZ(), collision.multFT0M());
       } else { // same-sign
-        // make sure to check for autocorrelations - only possible in same-sign correlations
+        // make sure to check for autocorrelations - only possible in same-sign correlations (if PID is correct)
         if (posIdTrigg == posIdAssoc && negIdTrigg == negIdAssoc) {
           // LOGF(info, "same v0 in SS correlation! %d %d", v0dataTrigg.v0Id(), v0dataAssoc.v0Id());
           registry.fill(HIST("hAutoCorrelation"), 0);