PWGLF: Improved PID selection and DCA fitter

PWGLF/Tasks/Resonances/highmasslambdasvx.cxx

@@ -95,12 +95,12 @@ struct highmasslambdasvx {
   Configurable<bool> additionalEvSel{"additionalEvSel", true, "additionalEvSel"};
   Configurable<bool> additionalEvSel2{"additionalEvSel2", false, "additionalEvSel2"};
   // events
+  Configurable<bool> cnfabsdca{"cnfabsdca", false, "Use Abs DCA for secondary vertex fitting"};
   Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"};
   Configurable<float> cfgCutCentralityMax{"cfgCutCentralityMax", 50.0f, "Accepted maximum Centrality"};
   Configurable<float> cfgCutCentralityMin{"cfgCutCentralityMin", 30.0f, "Accepted minimum Centrality"};
   // proton track cut
   Configurable<bool> ispTdifferentialDCA{"ispTdifferentialDCA", true, "is pT differential DCA"};
-  Configurable<bool> rejectPID{"rejectPID", true, "Reject PID"};
   Configurable<float> confMinRot{"confMinRot", 5.0 * TMath::Pi() / 6.0, "Minimum of rotation"};
   Configurable<float> confMaxRot{"confMaxRot", 7.0 * TMath::Pi() / 6.0, "Maximum of rotation"};
   Configurable<float> confRapidity{"confRapidity", 0.8, "cut on Rapidity"};
@@ -116,15 +116,15 @@ struct highmasslambdasvx {
   Configurable<int> cfgTPCcluster{"cfgTPCcluster", 70, "Number of TPC cluster"};
   Configurable<bool> ispTdepPID{"ispTdepPID", true, "pT dependent PID"};
   Configurable<float> nsigmaCutTPC{"nsigmacutTPC", 3.0, "Value of the TPC Nsigma cut"};
-  Configurable<float> nsigmaCutTOF{"nsigmacutTOF", 3.0, "Value of the TOF Nsigma cut"};
+  Configurable<float> nsigmaCutCombined{"nsigmaCutCombined", 2.0, "Combined nsigma cut"};
+  Configurable<float> nsigmaCutCombinedKaon{"nsigmaCutCombinedKaon", 1.5, "Combined nsigma cut Kaon"};
+  Configurable<float> nsigmaCutCombinedPion{"nsigmaCutCombinedPion", 2.0, "Combined nsigma cut Pion"};
   Configurable<float> nsigmaCutTPCPre{"nsigmacutTPCPre", 5.0, "Value of the TPC Nsigma cut Pre filter"};
   // Configs for V0
-  Configurable<float> ConfV0decaylengthmin{"ConfV0decaylengthmin", 2, "Minimum V0 decay length"};
   Configurable<float> ConfV0decaylengthmax{"ConfV0decaylengthmax", 100, "Maximum V0 decay length"};
   Configurable<float> ConfV0PtMin{"ConfV0PtMin", 0.f, "Minimum transverse momentum of V0"};
   Configurable<double> ConfV0DCADaughMax{"ConfV0DCADaughMax", 0.2f, "Maximum DCA between the V0 daughters"};
   Configurable<double> ConfV0CPAMin{"ConfV0CPAMin", 0.9998f, "Minimum CPA of V0"};
-  Configurable<float> ConfV0TranRadV0Min{"ConfV0TranRadV0Min", 1.5f, "Minimum transverse radius"};
   Configurable<float> ConfV0TranRadV0Max{"ConfV0TranRadV0Max", 100.f, "Maximum transverse radius"};
   Configurable<double> cMaxV0DCA{"cMaxV0DCA", 0.1, "Maximum V0 DCA to PV"};
   Configurable<float> cMaxV0LifeTime{"cMaxV0LifeTime", 20, "Maximum V0 life time"};
@@ -135,6 +135,10 @@ struct highmasslambdasvx {
   Configurable<float> ConfDaughTPCnclsMin{"ConfDaughTPCnclsMin", 50.f, "V0 Daugh sel: Min. nCls TPC"};
   Configurable<double> ConfDaughDCAMin{"ConfDaughDCAMin", 0.08f, "V0 Daugh sel:  Max. DCA Daugh to PV (cm)"};
   Configurable<float> ConfDaughPIDCuts{"ConfDaughPIDCuts", 3, "PID selections for KS0 daughters"};
+  // SVX cut
+  Configurable<bool> useCPAlambdac{"useCPAlambdac", true, "use CPA lambdac"};
+  Configurable<float> cutCPAlambdac{"cutCPAlambdac", 0.0, "cut CPA lambdac"};
+  Configurable<float> cutmaxctaulambdac{"cutmaxctaulambdac", 0.06, "cut max ctau lambdac"};
   // Fill strategy
   // Configurable<int> cfgSelectDaughterTopology{"cfgSelectDaughterTopology", 2, "Select daughter for topology"};
   // Mixed event
@@ -148,7 +152,7 @@ struct highmasslambdasvx {
   ConfigurableAxis configThnAxisSA{"configThnAxisSA", {100, -1, 1}, "SA"};
   ConfigurableAxis configThnAxisPhiminusPsi{"configThnAxisPhiminusPsi", {6, 0.0, TMath::Pi()}, "#phi - #psi"};
   ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {1, 30., 50}, "Centrality"};
-  ConfigurableAxis configThnAxisDecayLength{"configThnAxisDecayLength", {100, 0.0, 0.1}, "Decay length"};
+  ConfigurableAxis configThnAxisDecayLength{"configThnAxisDecayLength", {60, 0.0, 0.06}, "Decay length"};
 
   Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
   Filter centralityFilter = (nabs(aod::cent::centFT0C) < cfgCutCentralityMax && nabs(aod::cent::centFT0C) > cfgCutCentralityMin);
@@ -162,7 +166,7 @@ struct highmasslambdasvx {
 
   using TrackCandidates = soa::Filtered<soa::Join<aod::TracksIU, aod::TracksCovIU, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullPr, aod::pidTOFFullPr, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullKa, aod::pidTOFFullKa>>;
   using AllTrackCandidates = soa::Join<aod::TracksIU, aod::TracksCovIU, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullPi>;
-  using ResoV0s = aod::V0Datas;
+  using ResoV0s = soa::Join<aod::V0Datas, aod::V0Covs, aod::V0DauCovs>;
 
   SliceCache cache;
   // Partition<TrackCandidates> posTracks = aod::track::signed1Pt > cfgCutCharge;
@@ -199,11 +203,10 @@ struct highmasslambdasvx {
     histos.add("hSparseV2SASameEventRotational_V2_IOP", "hSparseV2SASameEventRotational_V2_IOP", HistType::kTHnSparseF, {thnAxisInvMass, ptAxis, thnAxisPhiminusPsi, dcaAxis, decaylengthAxis, occupancyAxis});
     histos.add("hSparseV2SAMixedEvent_V2_IOP", "hSparseV2SAMixedEvent_V2_IOP", HistType::kTHnSparseF, {thnAxisInvMass, ptAxis, thnAxisPhiminusPsi, dcaAxis, decaylengthAxis, occupancyAxis});
 
-    histos.add("hRejectPID", "hRejectPID", kTH1F, {{2, 0.0f, 2.0f}});
     histos.add("hV0decaylength", "hV0decaylength", kTH1F, {{1000, 0.0f, 1000.0f}});
     histos.add("hMomCorr", "hMomCorr", kTH3F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}, {8, 0.0f, 80.0f}});
     histos.add("hInvMassKs0", "hInvMassKs0", kTH1F, {{200, 0.4f, 0.6f}});
-    histos.add("hchi2PCA", "hchi2PCA", kTH1F, {{10000, 0.0f, 10000.f}});
+    histos.add("hchi2PCA", "hchi2PCA", kTH1F, {{1000, 0.0f, 1.f}});
     histos.add("hFTOCvsTPCNoCut", "Mult correlation FT0C vs. TPC without any cut", kTH2F, {{80, 0.0f, 80.0f}, {100, -0.5f, 5999.5f}});
     histos.add("hFTOCvsTPC", "Mult correlation FT0C vs. TPC", kTH2F, {{80, 0.0f, 80.0f}, {100, -0.5f, 5999.5f}});
     histos.add("hFTOCvsTPCSelected", "Mult correlation FT0C vs. TPC after selection", kTH2F, {{80, 0.0f, 80.0f}, {100, -0.5f, 5999.5f}});
@@ -221,6 +224,7 @@ struct highmasslambdasvx {
     histos.add("hDcaxy", "Dcaxy distribution", kTH1F, {{1000, -0.5f, 0.5f}});
     histos.add("hsignprotonDca", "hsignprotonDca", kTH1F, {{2000, -1.0f, 1.0f}});
     histos.add("hunsignprotonDca", "hunsignprotonDca", kTH1F, {{2000, -1.0f, 1.0f}});
+    histos.add("hsignCPALambdac", "hsignCPALambdac", kTH1F, {{2000, -1.0f, 1.0f}});
 
     histos.add("hPsiFT0C", "PsiFT0C", kTH3F, {centAxis, phiAxis, occupancyAxis});
     histos.add("hPsiFT0A", "PsiFT0A", kTH3F, {centAxis, phiAxis, occupancyAxis});
@@ -241,7 +245,7 @@ struct highmasslambdasvx {
     df.setMaxDZIni(4);
     df.setMinParamChange(1.e-3);
     df.setMinRelChi2Change(0.9);
-    df.setUseAbsDCA(true);
+    df.setUseAbsDCA(cnfabsdca);
     df.setWeightedFinalPCA(true);
     ccdb->setURL(ccdburl);
     ccdb->setCaching(true);
@@ -268,108 +272,48 @@ struct highmasslambdasvx {
   }
 
   template <typename T>
-  bool selectionPIDpTdependent(const T& candidate)
+  bool selectionPID(const T& candidate)
   {
-    if (candidate.p() <= 0.5 && TMath::Abs(candidate.tpcNSigmaPr()) < 5.0) {
+    if (candidate.pt() < 0.7 && !candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC) {
       return true;
     }
-    if (candidate.p() > 0.5 && candidate.p() <= 0.8 && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0) {
+    if (candidate.pt() >= 0.7 && !candidate.hasTOF() && candidate.pt() < 0.8 && candidate.tpcNSigmaPr() > -2.0 && candidate.tpcNSigmaPr() < nsigmaCutTPC) {
       return true;
     }
-    if (candidate.p() > 0.8 && candidate.hasTOF() && TMath::Sqrt(candidate.tofNSigmaPr() * candidate.tofNSigmaPr() + candidate.tpcNSigmaPr() * candidate.tpcNSigmaPr()) < nsigmaCutTOF) {
+    if (candidate.pt() < 0.8 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC && TMath::Abs(candidate.tofNSigmaPr()) < 10.0) {
       return true;
     }
-    return false;
-  }
-
-  template <typename T>
-  bool selectionPID(const T& candidate)
-  {
-    if (candidate.hasTOF()) {
-      if (candidate.pt() < 0.8 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && TMath::Abs(candidate.tofNSigmaPr()) < 10.0) {
-        return true;
-      }
-      if (candidate.pt() >= 0.8 && candidate.pt() < 3.0 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && TMath::Abs(candidate.tofNSigmaPr()) < 3.0) {
-        return true;
-      }
-      if (candidate.pt() >= 3.0 && candidate.pt() < 4.0 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && candidate.tofNSigmaPr() > -2.0 && candidate.tofNSigmaPr() < 3.0) {
-        return true;
-      }
-      if (candidate.pt() >= 4.0 && candidate.pt() < 5.0 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && candidate.tofNSigmaPr() > -1.5 && candidate.tofNSigmaPr() < 3.0) {
-        return true;
-      }
-      if (candidate.pt() >= 5.0 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && candidate.tofNSigmaPr() > -1.5 && candidate.tofNSigmaPr() < 3.0) {
-        return true;
-      }
+    if (candidate.pt() >= 0.8 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr() * candidate.tpcNSigmaPr() + candidate.tofNSigmaPr() * candidate.tofNSigmaPr()) < 2.0 * nsigmaCutCombined * nsigmaCutCombined && TMath::Abs(candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa() + candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) > 2.0 * nsigmaCutCombinedKaon * nsigmaCutCombinedKaon && TMath::Abs(candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) > 2.0 * nsigmaCutCombinedPion * nsigmaCutCombinedPion) {
+      return true;
     }
-    if (!candidate.hasTOF()) {
-      if (candidate.pt() < 0.7 && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0) {
-        return true;
-      }
-      if (candidate.pt() >= 0.7 && candidate.pt() < 0.8 && candidate.tpcNSigmaPr() > -2.0 && candidate.tpcNSigmaPr() < 3.0) {
-        return true;
-      }
-      if (candidate.pt() >= 0.8 && candidate.tpcNSigmaPr() > -1.5 && candidate.tpcNSigmaPr() < 3.0) {
-        return true;
-      }
+    if (candidate.pt() >= 0.8 && !candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutCombined && TMath::Abs(candidate.tpcNSigmaKa()) > nsigmaCutCombinedKaon && TMath::Abs(candidate.tpcNSigmaPi()) > nsigmaCutCombinedPion) {
+      return true;
     }
+
     return false;
   }
 
   template <typename T>
   bool selectionPIDNew(const T& candidate)
   {
-    if (candidate.pt() < 0.7 && !candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0) {
-      return true;
-    }
-    if (candidate.pt() >= 0.7 && !candidate.hasTOF() && candidate.pt() < 0.8 && candidate.tpcNSigmaPr() > -2.0 && candidate.tpcNSigmaPr() < 3.0) {
-      return true;
-    }
-    if (candidate.pt() < 0.8 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && TMath::Abs(candidate.tofNSigmaPr()) < 10.0) {
-      return true;
-    }
-    if (candidate.pt() >= 0.8 && candidate.pt() < 1.0 && !candidate.hasTOF() && candidate.tpcNSigmaPr() > -1.5 && candidate.tpcNSigmaPr() < 3.0) {
-      return true;
-    }
-    if (candidate.pt() >= 0.8 && candidate.pt() < 3.0 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && TMath::Abs(candidate.tofNSigmaPr()) < 3.0) {
+    if (candidate.pt() < 0.7 && !candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC) {
       return true;
     }
-    if (candidate.pt() >= 3.0 && candidate.pt() < 4.0 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < 3.0 && candidate.tofNSigmaPr() > -2.0 && candidate.tofNSigmaPr() < 4.0) {
+    if (candidate.pt() >= 0.7 && !candidate.hasTOF() && candidate.pt() < 0.8 && candidate.tpcNSigmaPr() > -2.0 && candidate.tpcNSigmaPr() < nsigmaCutTPC) {
       return true;
     }
-    if (candidate.pt() >= 4.0 && candidate.pt() < 5.0 && candidate.hasTOF() && candidate.tofNSigmaPr() > -2.0) {
+    if (candidate.pt() < 0.8 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC && TMath::Abs(candidate.tofNSigmaPr()) < 10.0) {
       return true;
     }
-    if (candidate.pt() >= 5.0 && candidate.hasTOF() && candidate.tofNSigmaPr() > -1.5) {
+    if (candidate.pt() >= 0.8 && candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaPr() * candidate.tpcNSigmaPr() + candidate.tofNSigmaPr() * candidate.tofNSigmaPr()) < 2.0 * nsigmaCutCombined * nsigmaCutCombined && TMath::Abs(candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa() + candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) > 2.0 * nsigmaCutCombinedKaon * nsigmaCutCombinedKaon && TMath::Abs(candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi()) > 2.0 * nsigmaCutCombinedPion * nsigmaCutCombinedPion) {
       return true;
     }
     return false;
   }
 
-  template <typename T>
-  bool RejectPion(const T& candidate)
-  {
-    if (candidate.p() > 1.0 && candidate.p() < 2.0 && !candidate.hasTOF() && candidate.tpcNSigmaPi() < 2) {
-      return false;
-    }
-    return true;
-  }
-
-  template <typename T>
-  bool RejectKaon(const T& candidate)
-  {
-    if (candidate.p() > 1.0 && candidate.p() < 2.0 && !candidate.hasTOF() && candidate.tpcNSigmaKa() < 2) {
-      return false;
-    }
-    return true;
-  }
-
   template <typename Collision, typename V0>
   bool SelectionV0(Collision const& collision, V0 const& candidate)
   {
-    if (fabs(candidate.dcav0topv()) > cMaxV0DCA) {
-      return false;
-    }
     const float pT = candidate.pt();
     const std::vector<float> decVtx = {candidate.x(), candidate.y(), candidate.z()};
     const float tranRad = candidate.v0radius();
@@ -383,22 +327,22 @@ struct highmasslambdasvx {
     if (pT < ConfV0PtMin) {
       return false;
     }
-    if ((decaylength < ConfV0decaylengthmin) || (decaylength > ConfV0decaylengthmax)) {
+    if (fabs(CtauK0s) < 2.0 || fabs(CtauK0s) > cMaxV0LifeTime || candidate.mK0Short() < lowmasscutks0 || candidate.mK0Short() > highmasscutks0) {
       return false;
     }
     if (dcaDaughv0 > ConfV0DCADaughMax) {
       return false;
     }
-    if (cpav0 < ConfV0CPAMin) {
+    if (fabs(candidate.dcav0topv()) > cMaxV0DCA) {
       return false;
     }
-    if (tranRad < ConfV0TranRadV0Min) {
+    if (cpav0 < ConfV0CPAMin) {
       return false;
     }
-    if (tranRad > ConfV0TranRadV0Max) {
+    if (decaylength > ConfV0decaylengthmax) {
       return false;
     }
-    if (fabs(CtauK0s) > cMaxV0LifeTime || candidate.mK0Short() < lowmasscutks0 || candidate.mK0Short() > highmasscutks0) {
+    if (tranRad > ConfV0TranRadV0Max) {
       return false;
     }
     return true;
@@ -523,14 +467,7 @@ struct highmasslambdasvx {
       if (!ispTdepPID && !selectionPID(track1)) {
         continue;
       }
-      if (!ispTdepPID && rejectPID && !RejectPion(track1)) {
-        histos.fill(HIST("hRejectPID"), 0.5);
-        continue;
-      }
-      if (!ispTdepPID && rejectPID && !RejectKaon(track1)) {
-        histos.fill(HIST("hRejectPID"), 1.5);
-        continue;
-      }
+
       histos.fill(HIST("hMomCorr"), track1.p() / track1.sign(), track1.p() - track1.tpcInnerParam(), centrality);
       histos.fill(HIST("hEta"), track1.eta());
       histos.fill(HIST("hDcaz"), track1.dcaZ());
@@ -580,11 +517,16 @@ struct highmasslambdasvx {
         v0y = v0.y();
         v0z = v0.z();
         const std::array<float, 3> vertexV0 = {v0x, v0y, v0z};
+
         v0px = v0.px();
         v0py = v0.py();
         v0pz = v0.pz();
         const std::array<float, 3> momentumV0 = {v0px, v0py, v0pz};
 
+        std::array<float, 6> covV0Pos = {0.};
+        for (int i = 0; i < 6; i++) {
+          covV0Pos[i] = v0.positionCovMat()[i];
+        }
         // dcaV0dau = v0.dcaV0daughters();
         // dcaPosToPV = v0.dcapostopv();
         // dcaNegToPV = v0.dcanegtopv();
@@ -594,7 +536,8 @@ struct highmasslambdasvx {
         trackParCovV0DaughNeg.propagateTo(negTrackX, bz); // propagate the track to the X closest to the V0 vertex
 
         // we build the neutral track to then build the cascade
-        auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, {0, 0, 0, 0, 0, 0}, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
+        // auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, {0, 0, 0, 0, 0, 0}, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
+        auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, covV0Pos, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
         std::array<float, 3> pVecV0 = {0., 0., 0.};
         std::array<float, 3> pVecBach = {0., 0., 0.};
         std::array<float, 3> pVecCand = {0., 0., 0.};
@@ -620,6 +563,9 @@ struct highmasslambdasvx {
 
         const auto& secondaryVertex = df.getPCACandidate();
         auto chi2PCA = df.getChi2AtPCACandidate();
+        if (chi2PCA > 1) {
+          continue;
+        }
         // auto covMatrixPCA = df.calcPCACovMatrixFlat();
         if (firstprimarytrack == 0) {
           histos.fill(HIST("hchi2PCA"), chi2PCA);
@@ -647,21 +593,28 @@ struct highmasslambdasvx {
         }
         double protonimpactparameter = impactParameter1.getY();
         // double kshortimpactparameter=impactParameter0.getY();
-        if (firstprimarytrack == 0) {
-          histos.fill(HIST("hDcaxy"), track1.dcaXY());
-          histos.fill(HIST("hunsignprotonDca"), protonimpactparameter);
-        }
-        double decaylengthx = collision.posX() - secondaryVertex[0];
-        double decaylengthy = collision.posY() - secondaryVertex[1];
-        double decaylengthz = collision.posZ() - secondaryVertex[2];
+
+        double decaylengthx = secondaryVertex[0] - collision.posX();
+        double decaylengthy = secondaryVertex[1] - collision.posY();
+        double decaylengthz = secondaryVertex[2] - collision.posZ();
         double decaylength = TMath::Sqrt(decaylengthx * decaylengthx + decaylengthy * decaylengthy + decaylengthz * decaylengthz);
-        double lambdaclifetime = decaylength;
-        double anglesign = decaylengthx * Lambdac.X() + decaylengthy * Lambdac.Y() + decaylengthz * Lambdac.Z();
+        double lambdaclifetime = (decaylength * Lambdac.M()) / Lambdac.P();
+        if (lambdaclifetime > cutmaxctaulambdac) {
+          continue;
+        }
+        double anglesign = decaylengthx * Lambdac.Px() + decaylengthy * Lambdac.Py() + decaylengthz * Lambdac.Pz();
+        double CPAlambdac = anglesign / (decaylength * Lambdac.P());
         anglesign = anglesign / TMath::Abs(anglesign);
         auto signprotonimpactparameter = protonimpactparameter * anglesign;
         if (firstprimarytrack == 0) {
+          histos.fill(HIST("hDcaxy"), track1.dcaXY());
+          histos.fill(HIST("hunsignprotonDca"), protonimpactparameter);
+          histos.fill(HIST("hsignCPALambdac"), TMath::Abs(CPAlambdac));
           histos.fill(HIST("hsignprotonDca"), signprotonimpactparameter);
         }
+        if (useCPAlambdac && TMath::Abs(CPAlambdac) < cutCPAlambdac) {
+          continue;
+        }
         firstprimarytrack = firstprimarytrack + 1;
         auto phiminuspsi = GetPhiInRange(Lambdac.Phi() - psiFT0C);
         v2 = TMath::Cos(2.0 * phiminuspsi);