PWGHF: Add dimuon channel to X task (#78)

lucamicheletti93 · Luca Micheletti · vkucera · web-flow · commit 983f119a62ab · 2021-09-27T20:03:49.000+02:00
* Adding muon channel to X

* Running dielectron and dimuon channels separately

* Implementing comments

* Small fix

* Small Fix

* Implementing comments

* Fixing the issue with the X pdg code

* Fixing the issue with the X pdg code

* Modification to select candidates in both the channels

* Checking the Jpsi is in the decay chain

* Checking the Jpsi is in the decay chain

* Small fix

* Moving X functions from hf_cand_prong3 to hf_cand_x

* Small fix

* Fix selector

Co-authored-by: Luca Micheletti &lt;lmichele@alicecerno2.cern.ch&gt;
Co-authored-by: Vít Kučera &lt;vit.kucera@cern.ch&gt;
diff --git a/PWGHF/DataModel/HFCandidateSelectionTables.h b/PWGHF/DataModel/HFCandidateSelectionTables.h
@@ -85,10 +85,11 @@ DECLARE_SOA_TABLE(HFSelBPlusToD0PiCandidate, "AOD", "HFSELBPLUSCAND", //!
 
 namespace hf_selcandidate_x
 {
-DECLARE_SOA_COLUMN(IsSelXToJpsiPiPi, isSelXToJpsiPiPi, int); //!
+DECLARE_SOA_COLUMN(IsSelXToJpsiToEEPiPi, isSelXToJpsiToEEPiPi, int);     //!
+DECLARE_SOA_COLUMN(IsSelXToJpsiToMuMuPiPi, isSelXToJpsiToMuMuPiPi, int); //!
 } // namespace hf_selcandidate_x
 DECLARE_SOA_TABLE(HFSelXToJpsiPiPiCandidate, "AOD", "HFSELXCAND", //!
-                  hf_selcandidate_x::IsSelXToJpsiPiPi);
+                  hf_selcandidate_x::IsSelXToJpsiToEEPiPi, hf_selcandidate_x::IsSelXToJpsiToMuMuPiPi);
 } // namespace o2::aod
 
 namespace o2::aod
diff --git a/PWGHF/DataModel/HFSecondaryVertex.h b/PWGHF/DataModel/HFSecondaryVertex.h
@@ -702,32 +702,6 @@ auto InvMassXicToPiKP(const T& candidate)
   return InvMassLcpiKp(candidate);
 }
 
-// X → Jpsi π+ π-
-// TODO: add pdg code for X (9920443), temporarily hardcode mass here:
-float massX = 3.872; // replace this with: "RecoDecay::getMassPDG(9920443)" when pdg is added
-template <typename T>
-auto CtX(const T& candidate)
-{
-  return candidate.ct(massX);
-}
-
-template <typename T>
-auto YX(const T& candidate)
-{
-  return candidate.y(massX);
-}
-
-template <typename T>
-auto EX(const T& candidate)
-{
-  return candidate.e(massX);
-}
-template <typename T>
-auto InvMassXToJpsiPiPi(const T& candidate)
-{
-  return candidate.m(array{RecoDecay::getMassPDG(443), RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kPiPlus)});
-}
-
 } // namespace hf_cand_prong3
 
 // 3-prong decay candidate table
@@ -800,9 +774,36 @@ DECLARE_SOA_COLUMN(OriginMCGen, originMCGen, int8_t);               // particle
 DECLARE_SOA_COLUMN(FlagMCDecayChanRec, flagMCDecayChanRec, int8_t); // resonant decay channel flag, reconstruction level
 DECLARE_SOA_COLUMN(FlagMCDecayChanGen, flagMCDecayChanGen, int8_t); // resonant decay channel flag, generator level
 // mapping of decay types
-enum DecayType { XToJpsiPiPi = 0 }; // move this to a dedicated cascade namespace in the future?
+enum DecayType { XToJpsiToEEPiPi = 0,
+                 XToJpsiToMuMuPiPi }; // move this to a dedicated cascade namespace in the future?
 } // namespace hf_cand_x
 
+// X → Jpsi π+ π-
+// TODO: add pdg code for X (9920443), temporarily hardcode mass here:
+float massX = 3.872; // replace this with: "RecoDecay::getMassPDG(9920443)" when pdg is added
+template <typename T>
+auto CtX(const T& candidate)
+{
+  return candidate.ct(massX);
+}
+
+template <typename T>
+auto YX(const T& candidate)
+{
+  return candidate.y(massX);
+}
+
+template <typename T>
+auto EX(const T& candidate)
+{
+  return candidate.e(massX);
+}
+template <typename T>
+auto InvMassXToJpsiPiPi(const T& candidate)
+{
+  return candidate.m(array{RecoDecay::getMassPDG(443), RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kPiPlus)});
+}
+
 // declare dedicated X candidate table
 DECLARE_SOA_TABLE(HfCandXBase, "AOD", "HFCANDXBASE",
                   // general columns
diff --git a/PWGHF/TableProducer/HFCandidateCreatorX.cxx b/PWGHF/TableProducer/HFCandidateCreatorX.cxx
@@ -14,6 +14,7 @@
 /// \note Adapted from HFCandidateCreator3Prong
 ///
 /// \author Rik Spijkers <r.spijkers@students.uu.nl>, Utrecht University
+/// \author Luca Micheletti <luca.micheletti@to.infn.it>, INFN
 
 #include "Framework/AnalysisTask.h"
 #include "DetectorsVertexing/DCAFitterN.h"
@@ -24,10 +25,10 @@
 #include "PWGHF/DataModel/HFCandidateSelectionTables.h"
 
 using namespace o2;
+using namespace o2::aod;
 using namespace o2::framework;
 using namespace o2::aod::hf_cand;
 using namespace o2::aod::hf_cand_prong2;
-using namespace o2::aod::hf_cand_prong3;
 using namespace o2::aod::hf_cand_x;
 using namespace o2::framework::expressions;
 
@@ -52,10 +53,13 @@ struct HFCandidateCreatorX {
   Configurable<double> ptPionMin{"ptPionMin", 0.15, "minimum pion pT threshold (GeV/c)"};
   Configurable<bool> b_dovalplots{"b_dovalplots", true, "do validation plots"};
 
-  OutputObj<TH1F> hMassJpsi{TH1F("hMassJpsi", "J/#psi candidates;inv. mass (#e+ e-) (GeV/#it{c}^{2});entries", 500, 0., 5.)};
+  OutputObj<TH1F> hMassJpsiToEE{TH1F("hMassJpsiToEE", "J/#psi candidates;inv. mass (e^{#plus} e^{#minus}) (GeV/#it{c}^{2});entries", 500, 0., 5.)};
+  OutputObj<TH1F> hMassJpsiToMuMu{TH1F("hMassJpsiToMuMu", "J/#psi candidates;inv. mass (#mu^{#plus} #mu^{#minus}) (GeV/#it{c}^{2});entries", 500, 0., 5.)};
   OutputObj<TH1F> hPtJpsi{TH1F("hPtJpsi", "J/#psi candidates;candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)};
+  OutputObj<TH1F> hPtPion{TH1F("hPtPion", "#pi candidates;candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)};
   OutputObj<TH1F> hCPAJpsi{TH1F("hCPAJpsi", "J/#psi candidates;cosine of pointing angle;entries", 110, -1.1, 1.1)};
-  OutputObj<TH1F> hMassX{TH1F("hMassX", "3-prong candidates;inv. mass (J/#psi #pi+ #pi-) (GeV/#it{c}^{2});entries", 500, 0., 5.)};
+  OutputObj<TH1F> hMassXToJpsiToEEPiPi{TH1F("hMassXToJpsiToEEPiPi", "3-prong candidates;inv. mass (J/#psi (#rightarrow e+ e-) #pi+ #pi-) (GeV/#it{c}^{2});entries", 500, 0., 5.)};
+  OutputObj<TH1F> hMassXToJpsiToMuMuPiPi{TH1F("hMassXToJpsiToMuMuPiPi", "3-prong candidates;inv. mass (J/#psi (#rightarrow #mu+ #mu-) #pi+ #pi-) (GeV/#it{c}^{2});entries", 500, 0., 5.)};
   OutputObj<TH1F> hCovPVXX{TH1F("hCovPVXX", "3-prong candidates;XX element of cov. matrix of prim. vtx position (cm^{2});entries", 100, 0., 1.e-4)};
   OutputObj<TH1F> hCovSVXX{TH1F("hCovSVXX", "3-prong candidates;XX element of cov. matrix of sec. vtx position (cm^{2});entries", 100, 0., 0.2)};
 
@@ -65,7 +69,7 @@ struct HFCandidateCreatorX {
 
   Configurable<int> d_selectionFlagJpsi{"d_selectionFlagJpsi", 1, "Selection Flag for Jpsi"};
   Configurable<double> cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"};
-  Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEE >= d_selectionFlagJpsi);
+  Filter filterSelectCandidates = (aod::hf_selcandidate_jpsi::isSelJpsiToEE >= d_selectionFlagJpsi || aod::hf_selcandidate_jpsi::isSelJpsiToMuMu >= d_selectionFlagJpsi);
 
   void process(aod::Collision const& collision,
                soa::Filtered<soa::Join<
@@ -95,13 +99,18 @@ struct HFCandidateCreatorX {
 
     // loop over Jpsi candidates
     for (auto& jpsiCand : jpsiCands) {
-      if (!(jpsiCand.hfflag() & 1 << JpsiToEE)) {
+      if (!(jpsiCand.hfflag() & 1 << hf_cand_prong2::DecayType::JpsiToEE) && !(jpsiCand.hfflag() & 1 << hf_cand_prong2::DecayType::JpsiToMuMu)) {
         continue;
       }
       if (cutYCandMax >= 0. && std::abs(YJpsi(jpsiCand)) > cutYCandMax) {
         continue;
       }
-      hMassJpsi->Fill(InvMassJpsiToEE(jpsiCand));
+      if (jpsiCand.isSelJpsiToEE() > 0) {
+        hMassJpsiToEE->Fill(InvMassJpsiToEE(jpsiCand));
+      }
+      if (jpsiCand.isSelJpsiToMuMu() > 0) {
+        hMassJpsiToMuMu->Fill(InvMassJpsiToMuMu(jpsiCand));
+      }
       hPtJpsi->Fill(jpsiCand.pt());
       hCPAJpsi->Fill(jpsiCand.cpa());
       // create Jpsi track to pass to DCA fitter; use cand table + rebuild vertex
@@ -138,6 +147,7 @@ struct HFCandidateCreatorX {
         if (trackPos.pt() < ptPionMin) {
           continue;
         }
+        hPtPion->Fill(trackPos.pt());
 
         // loop over pi- candidates
         for (auto& trackNeg : tracks) {
@@ -150,6 +160,7 @@ struct HFCandidateCreatorX {
           if (trackNeg.pt() < ptPionMin) {
             continue;
           }
+          hPtPion->Fill(trackNeg.pt());
 
           auto trackParVarPos = getTrackParCov(trackPos);
           auto trackParVarNeg = getTrackParCov(trackNeg);
@@ -190,7 +201,13 @@ struct HFCandidateCreatorX {
           auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta));
           auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.));
 
-          int hfFlag = 1 << XToJpsiPiPi;
+          int hfFlag = 0;
+          if (TESTBIT(jpsiCand.hfflag(), hf_cand_prong2::DecayType::JpsiToMuMu)) {
+            SETBIT(hfFlag, hf_cand_x::DecayType::XToJpsiToMuMuPiPi); // dimuon channel
+          }
+          if (TESTBIT(jpsiCand.hfflag(), hf_cand_prong2::DecayType::JpsiToEE)) {
+            SETBIT(hfFlag, hf_cand_x::DecayType::XToJpsiToEEPiPi); // dielectron channel
+          }
 
           // fill the candidate table for the X here:
           rowCandidateBase(collision.globalIndex(),
@@ -209,7 +226,12 @@ struct HFCandidateCreatorX {
           // calculate invariant mass
           auto arrayMomenta = array{pvecJpsi, pvecPos, pvecNeg};
           massJpsiPiPi = RecoDecay::M(std::move(arrayMomenta), array{massJpsi, massPi, massPi});
-          hMassX->Fill(massJpsiPiPi);
+          if (jpsiCand.isSelJpsiToEE() > 0) {
+            hMassXToJpsiToEEPiPi->Fill(massJpsiPiPi);
+          }
+          if (jpsiCand.isSelJpsiToMuMu() > 0) {
+            hMassXToJpsiToMuMuPiPi->Fill(massJpsiPiPi);
+          }
         } // pi- loop
       }   // pi+ loop
     }     // Jpsi loop
@@ -233,6 +255,8 @@ struct HFCandidateCreatorXMC {
                aod::McParticles const& particlesMC)
   {
     int indexRec = -1;
+    int pdgCodeX = pdg::Code::kX3872;
+    int pdgCodeJpsi = pdg::Code::kJpsi;
     int8_t sign = 0;
     int8_t flag = 0;
     int8_t origin = 0;
@@ -245,16 +269,32 @@ struct HFCandidateCreatorXMC {
       origin = 0;
       channel = 0;
       auto jpsiTrack = candidate.index0();
+      auto daughterPosJpsi = jpsiTrack.index0_as<aod::BigTracksMC>();
+      auto daughterNegJpsi = jpsiTrack.index1_as<aod::BigTracksMC>();
+      auto arrayJpsiDaughters = array{daughterPosJpsi, daughterNegJpsi};
       auto arrayDaughters = array{candidate.index1_as<aod::BigTracksMC>(),
                                   candidate.index2_as<aod::BigTracksMC>(),
-                                  jpsiTrack.index0_as<aod::BigTracksMC>(),
-                                  jpsiTrack.index1_as<aod::BigTracksMC>()};
+                                  daughterPosJpsi,
+                                  daughterNegJpsi};
 
       // X → J/ψ π+ π-
       //Printf("Checking X → J/ψ π+ π-");
-      indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, 9920443, array{+kPiPlus, -kPiPlus, +kElectron, -kElectron}, true, &sign, 2);
+      indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayJpsiDaughters, pdg::Code::kJpsi, array{+kElectron, -kElectron}, true);
       if (indexRec > -1) {
-        flag = 1 << XToJpsiPiPi;
+        indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, pdgCodeX, array{+kPiPlus, -kPiPlus, +kElectron, -kElectron}, true, &sign, 2);
+        if (indexRec > -1) {
+          flag = 1 << hf_cand_x::DecayType::XToJpsiToEEPiPi;
+        }
+      }
+
+      if (flag == 0) {
+        indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayJpsiDaughters, pdg::Code::kJpsi, array{+kMuonPlus, -kMuonPlus}, true);
+        if (indexRec > -1) {
+          indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, pdgCodeX, array{+kPiPlus, -kPiPlus, +kMuonPlus, -kMuonPlus}, true, &sign, 2);
+          if (indexRec > -1) {
+            flag = 1 << hf_cand_x::DecayType::XToJpsiToMuMuPiPi;
+          }
+        }
       }
 
       // Check whether the particle is non-prompt (from a b quark).
@@ -275,13 +315,19 @@ struct HFCandidateCreatorXMC {
 
       // X → J/ψ π+ π-
       //Printf("Checking X → J/ψ π+ π-");
-      if (RecoDecay::isMatchedMCGen(particlesMC, particle, 9920443, array{443, +kPiPlus, -kPiPlus}, true)) {
+      if (RecoDecay::isMatchedMCGen(particlesMC, particle, pdgCodeX, array{pdgCodeJpsi, +kPiPlus, -kPiPlus}, true)) {
         // Match J/psi --> e+e-
         std::vector<int> arrDaughter;
-        RecoDecay::getDaughters(particlesMC, particle, &arrDaughter, array{443}, 1);
+        RecoDecay::getDaughters(particlesMC, particle, &arrDaughter, array{pdgCodeJpsi}, 1);
         auto jpsiCandMC = particlesMC.iteratorAt(arrDaughter[0]);
-        if (RecoDecay::isMatchedMCGen(particlesMC, jpsiCandMC, 443, array{+kElectron, -kElectron}, true)) {
-          flag = 1 << XToJpsiPiPi;
+        if (RecoDecay::isMatchedMCGen(particlesMC, jpsiCandMC, pdgCodeJpsi, array{+kElectron, -kElectron}, true)) {
+          flag = 1 << hf_cand_x::DecayType::XToJpsiToEEPiPi;
+        }
+
+        if (flag == 0) {
+          if (RecoDecay::isMatchedMCGen(particlesMC, jpsiCandMC, pdgCodeJpsi, array{+kMuonPlus, -kMuonPlus}, true)) {
+            flag = 1 << hf_cand_x::DecayType::XToJpsiToMuMuPiPi;
+          }
         }
       }
 
diff --git a/PWGHF/TableProducer/HFTreeCreatorXToJpsiPiPi.cxx b/PWGHF/TableProducer/HFTreeCreatorXToJpsiPiPi.cxx
@@ -26,6 +26,7 @@
 #include "ReconstructionDataFormats/DCA.h"
 
 using namespace o2;
+using namespace o2::aod;
 using namespace o2::framework;
 using namespace o2::aod::hf_cand_x;
 
@@ -173,14 +174,15 @@ struct HfTreeCreatorXTojpsipipi {
         }
       };
 
-      fillTable(0, candidate.isSelXToJpsiPiPi(), InvMassXToJpsiPiPi(candidate), CtX(candidate), YX(candidate));
+      fillTable(0, candidate.isSelXToJpsiToEEPiPi(), InvMassXToJpsiPiPi(candidate), CtX(candidate), YX(candidate));
+      fillTable(1, candidate.isSelXToJpsiToMuMuPiPi(), InvMassXToJpsiPiPi(candidate), CtX(candidate), YX(candidate));
     }
 
     // Filling particle properties
     float massX = 3.872;
     rowCandidateFullParticles.reserve(particles.size());
     for (auto& particle : particles) {
-      if (std::abs(particle.flagMCMatchGen()) == 1 << DecayType::XToJpsiPiPi) {
+      if (std::abs(particle.flagMCMatchGen()) == 1 << DecayType::XToJpsiToEEPiPi || std::abs(particle.flagMCMatchGen()) == 1 << DecayType::XToJpsiToMuMuPiPi) {
         rowCandidateFullParticles(
           particle.mcCollision().bcId(),
           particle.pt(),
diff --git a/PWGHF/TableProducer/HFXToJpsiPiPiCandidateSelector.cxx b/PWGHF/TableProducer/HFXToJpsiPiPiCandidateSelector.cxx
@@ -13,13 +13,15 @@
 /// \brief X(3872) selection task.
 /// \note Adapted from HFJpsiCandidateSelector.cxx
 /// \author Rik Spijkers <r.spijkers@students.uu.nl>, Utrecht University
+/// \author Luca Micheletti <luca.micheletti@to.infn.it>, INFN
 
 #include "Framework/runDataProcessing.h"
 #include "Framework/AnalysisTask.h"
 #include "PWGHF/Core/HFSelectorCuts.h"
 #include "PWGHF/DataModel/HFSecondaryVertex.h"
 #include "PWGHF/DataModel/HFCandidateSelectionTables.h"
 using namespace o2;
+using namespace o2::aod;
 using namespace o2::framework;
 using namespace o2::aod::hf_cand_x;
 using namespace o2::analysis;
@@ -191,16 +193,26 @@ struct HFXToJpsiPiPiCandidateSelector {
       auto trackPos = hfCandX.index1_as<aod::BigTracksPID>(); //positive daughter
       auto trackNeg = hfCandX.index2_as<aod::BigTracksPID>(); //negative daughter
 
+      int selJpsiToEE = 1;
+      int selJpsiToMuMu = 1;
+
       // check if flagged as X --> Jpsi Pi Pi
-      if (!(hfCandX.hfflag() & 1 << XToJpsiPiPi)) {
-        hfSelXToJpsiPiPiCandidate(0);
-        // Printf("X candidate selection failed at hfflag check");
+      if (!(hfCandX.hfflag() & 1 << hf_cand_x::DecayType::XToJpsiToEEPiPi)) {
+        selJpsiToEE = 0;
+      }
+
+      if (!(hfCandX.hfflag() & 1 << hf_cand_x::DecayType::XToJpsiToMuMuPiPi)) {
+        selJpsiToMuMu = 0;
+      }
+
+      if (selJpsiToEE == 0 && selJpsiToMuMu == 0) {
+        hfSelXToJpsiPiPiCandidate(0, 0);
         continue;
       }
 
       // daughter track validity selection
       if (!daughterSelection(trackPos) || !daughterSelection(trackNeg)) {
-        hfSelXToJpsiPiPiCandidate(0);
+        hfSelXToJpsiPiPiCandidate(0, 0);
         // Printf("X candidate selection failed at daughter selection");
         continue;
       }
@@ -209,18 +221,18 @@ struct HFXToJpsiPiPiCandidateSelector {
       //need to add special cuts (additional cuts on decay length and d0 norm)
 
       if (!selectionTopol(hfCandX, candJpsi, trackPos, trackNeg)) {
-        hfSelXToJpsiPiPiCandidate(0);
+        hfSelXToJpsiPiPiCandidate(0, 0);
         // Printf("X candidate selection failed at selection topology");
         continue;
       }
 
       if (selectionPID(trackPos) == 0 || selectionPID(trackNeg) == 0) {
-        hfSelXToJpsiPiPiCandidate(0);
+        hfSelXToJpsiPiPiCandidate(0, 0);
         // Printf("X candidate selection failed at selection PID");
         continue;
       }
 
-      hfSelXToJpsiPiPiCandidate(1);
+      hfSelXToJpsiPiPiCandidate(selJpsiToEE, selJpsiToMuMu);
       // Printf("X candidate selection successful, candidate should be selected");
     }
   }
diff --git a/PWGHF/Tasks/taskX.cxx b/PWGHF/Tasks/taskX.cxx
