[PWGHF] Add the possibility to run MC Gen only for HF derived data

PWGHF/TableProducer/CMakeLists.txt

@@ -1,4 +1,4 @@
 # Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 # See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 # All rights not expressly granted are reserved.
 #
@@ -100,6 +100,11 @@
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(candidate-creator-mc-gen
+                    SOURCES candidateCreatorMcGen.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 # Candidate selectors
 
 o2physics_add_dpl_workflow(candidate-selector-b0-to-d-pi
@@ -244,12 +249,12 @@
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter O2Physics::EventFilteringUtils
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(tree-creator-omegac0-to-omega-ka
                     SOURCES treeCreatorOmegacToOmegaKa.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(tree-creator-omegac0-to-omega-pi
                     SOURCES treeCreatorOmegacToOmegaPi.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)

PWGHF/TableProducer/candidateCreator2Prong.cxx

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -17,7 +17,7 @@
 /// \author Pengzhong Lu <pengzhong.lu@cern.ch>, GSI Darmstadt, USTC
 
 #ifndef HomogeneousField
 #define HomogeneousField
 #endif
 
 #include <memory>
@@ -49,6 +49,7 @@
 #include "PWGHF/Utils/utilsEvSelHf.h"
 #include "PWGHF/Utils/utilsPid.h"
 #include "PWGHF/Utils/utilsTrkCandHf.h"
+#include "PWGHF/Utils/utilsMcGen.h"
 
 using namespace o2;
 using namespace o2::analysis;
@@ -371,7 +372,7 @@
         kfpVertex.SetCovarianceMatrix(rowTrackIndexProng2.pvRefitSigmaX2(), rowTrackIndexProng2.pvRefitSigmaXY(), rowTrackIndexProng2.pvRefitSigmaY2(), rowTrackIndexProng2.pvRefitSigmaXZ(), rowTrackIndexProng2.pvRefitSigmaYZ(), rowTrackIndexProng2.pvRefitSigmaZ2());
       }
       kfpVertex.GetCovarianceMatrix(covMatrixPV);
       KFParticle KFPV(kfpVertex);
       registry.fill(HIST("hCovPVXX"), covMatrixPV[0]);
       registry.fill(HIST("hCovPVYY"), covMatrixPV[2]);
       registry.fill(HIST("hCovPVXZ"), covMatrixPV[3]);
@@ -692,10 +693,10 @@
   using McCollisionsFT0Cs = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Cs>;
   using McCollisionsFT0Ms = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Ms>;
   using McCollisionsCentFT0Ms = soa::Join<aod::McCollisions, aod::McCentFT0Ms>;
   PresliceUnsorted<McCollisionsNoCents> colPerMcCollision = aod::mccollisionlabel::mcCollisionId;
   PresliceUnsorted<McCollisionsFT0Cs> colPerMcCollisionFT0C = aod::mccollisionlabel::mcCollisionId;
   PresliceUnsorted<McCollisionsFT0Ms> colPerMcCollisionFT0M = aod::mccollisionlabel::mcCollisionId;
   Preslice<aod::McParticles> mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId;
 
   using BCsInfo = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
   HistogramRegistry registry{"registry"};
@@ -839,47 +840,7 @@
         }
         continue;
       }
-
-      // Match generated particles.
-      for (const auto& particle : mcParticlesPerMcColl) {
-        flag = 0;
-        origin = 0;
-        std::vector<int> idxBhadMothers{};
-        // Reject particles from background events
-        if (particle.fromBackgroundEvent() && rejectBackground) {
-          rowMcMatchGen(flag, origin, -1);
-          continue;
-        }
-
-        // D0(bar) → π± K∓
-        if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign)) {
-          flag = sign * (1 << DecayType::D0ToPiK);
-        }
-
-        // J/ψ → e+ e−
-        if (flag == 0) {
-          if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kJPsi, std::array{+kElectron, -kElectron}, true)) {
-            flag = 1 << DecayType::JpsiToEE;
-          }
-        }
-
-        // J/ψ → μ+ μ−
-        if (flag == 0) {
-          if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kJPsi, std::array{+kMuonPlus, -kMuonPlus}, true)) {
-            flag = 1 << DecayType::JpsiToMuMu;
-          }
-        }
-
-        // Check whether the particle is non-prompt (from a b quark).
-        if (flag != 0) {
-          origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers);
-        }
-        if (origin == RecoDecay::OriginType::NonPrompt) {
-          rowMcMatchGen(flag, origin, idxBhadMothers[0]);
-        } else {
-          rowMcMatchGen(flag, origin, -1);
-        }
-      }
+      hf_mc_gen::fillMcMatchGen2Prong(mcParticlesPerMcColl, rowMcMatchGen, rejectBackground);
     }
   }
 
@@ -917,6 +878,6 @@
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   return WorkflowSpec{
     adaptAnalysisTask<HfCandidateCreator2Prong>(cfgc, TaskName{"hf-candidate-creator-2prong"}),
     adaptAnalysisTask<HfCandidateCreator2ProngExpressions>(cfgc, TaskName{"hf-candidate-creator-2prong-expressions"})};
 }

PWGHF/TableProducer/candidateCreator3Prong.cxx

@@ -48,6 +48,7 @@
 #include "PWGHF/Utils/utilsBfieldCCDB.h"
 #include "PWGHF/Utils/utilsEvSelHf.h"
 #include "PWGHF/Utils/utilsTrkCandHf.h"
+#include "PWGHF/Utils/utilsMcGen.h"
 
 using namespace o2;
 using namespace o2::analysis;
@@ -1048,96 +1049,7 @@ struct HfCandidateCreator3ProngExpressions {
         }
         continue;
       }
-
-      // Match generated particles.
-      for (const auto& particle : mcParticlesPerMcColl) {
-        flag = 0;
-        origin = 0;
-        channel = 0;
-        arrDaughIndex.clear();
-        std::vector<int> idxBhadMothers{};
-        // Reject particles from background events
-        if (particle.fromBackgroundEvent() && rejectBackground) {
-          rowMcMatchGen(flag, origin, channel, -1);
-          continue;
-        }
-
-        // D± → π± K∓ π±
-        if (createDplus) {
-          if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) {
-            flag = sign * (1 << DecayType::DplusToPiKPi);
-          }
-        }
-
-        // Ds± → K± K∓ π± and D± → K± K∓ π±
-        if (flag == 0 && createDs) {
-          bool isDplus = false;
-          if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDS, std::array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) {
-            // DecayType::DsToKKPi is used to flag both Ds± → K± K∓ π± and D± → K± K∓ π±
-            // TODO: move to different and explicit flags
-            flag = sign * (1 << DecayType::DsToKKPi);
-          } else if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kDPlus, std::array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 2)) {
-            // DecayType::DsToKKPi is used to flag both Ds± → K± K∓ π± and D± → K± K∓ π±
-            // TODO: move to different and explicit flags
-            flag = sign * (1 << DecayType::DsToKKPi);
-            isDplus = true;
-          }
-          if (flag != 0) {
-            RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1);
-            if (arrDaughIndex.size() == 2) {
-              for (auto jProng = 0u; jProng < arrDaughIndex.size(); ++jProng) {
-                auto daughJ = mcParticles.rawIteratorAt(arrDaughIndex[jProng]);
-                arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode());
-              }
-              if ((arrPDGDaugh[0] == arrPDGResonantDPhiPi[0] && arrPDGDaugh[1] == arrPDGResonantDPhiPi[1]) || (arrPDGDaugh[0] == arrPDGResonantDPhiPi[1] && arrPDGDaugh[1] == arrPDGResonantDPhiPi[0])) {
-                channel = isDplus ? DecayChannelDToKKPi::DplusToPhiPi : DecayChannelDToKKPi::DsToPhiPi;
-              } else if ((arrPDGDaugh[0] == arrPDGResonantDKstarK[0] && arrPDGDaugh[1] == arrPDGResonantDKstarK[1]) || (arrPDGDaugh[0] == arrPDGResonantDKstarK[1] && arrPDGDaugh[1] == arrPDGResonantDKstarK[0])) {
-                channel = isDplus ? DecayChannelDToKKPi::DplusToK0starK : DecayChannelDToKKPi::DsToK0starK;
-              }
-            }
-          }
-        }
-
-        // Λc± → p± K∓ π±
-        if (flag == 0 && createLc) {
-          if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kLambdaCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) {
-            flag = sign * (1 << DecayType::LcToPKPi);
-
-            // Flagging the different Λc± → p± K∓ π± decay channels
-            RecoDecay::getDaughters(particle, &arrDaughIndex, std::array{0}, 1);
-            if (arrDaughIndex.size() == 2) {
-              for (auto jProng = 0u; jProng < arrDaughIndex.size(); ++jProng) {
-                auto daughJ = mcParticles.rawIteratorAt(arrDaughIndex[jProng]);
-                arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode());
-              }
-              if ((arrPDGDaugh[0] == arrPDGResonant1[0] && arrPDGDaugh[1] == arrPDGResonant1[1]) || (arrPDGDaugh[0] == arrPDGResonant1[1] && arrPDGDaugh[1] == arrPDGResonant1[0])) {
-                channel = 1;
-              } else if ((arrPDGDaugh[0] == arrPDGResonant2[0] && arrPDGDaugh[1] == arrPDGResonant2[1]) || (arrPDGDaugh[0] == arrPDGResonant2[1] && arrPDGDaugh[1] == arrPDGResonant2[0])) {
-                channel = 2;
-              } else if ((arrPDGDaugh[0] == arrPDGResonant3[0] && arrPDGDaugh[1] == arrPDGResonant3[1]) || (arrPDGDaugh[0] == arrPDGResonant3[1] && arrPDGDaugh[1] == arrPDGResonant3[0])) {
-                channel = 3;
-              }
-            }
-          }
-        }
-
-        // Ξc± → p± K∓ π±
-        if (flag == 0 && createXic) {
-          if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kXiCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) {
-            flag = sign * (1 << DecayType::XicToPKPi);
-          }
-        }
-
-        // Check whether the particle is non-prompt (from a b quark).
-        if (flag != 0) {
-          origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers);
-        }
-        if (origin == RecoDecay::OriginType::NonPrompt) {
-          rowMcMatchGen(flag, origin, channel, idxBhadMothers[0]);
-        } else {
-          rowMcMatchGen(flag, origin, channel, -1);
-        }
-      }
+      hf_mc_gen::fillMcMatchGen3Prong(mcParticlesPerMcColl, rowMcMatchGen, rejectBackground, createDplus, createDs, createLc, createXic);
     }
   }
 

PWGHF/TableProducer/candidateCreatorB0.cxx

@@ -15,6 +15,10 @@
 ///
 /// \author Alexandre Bigot <alexandre.bigot@cern.ch>, IPHC Strasbourg
 
+#include <vector>
+#include <string>
+#include <memory>
+
 #include "CommonConstants/PhysicsConstants.h"
 #include "DCAFitter/DCAFitterN.h"
 #include "Framework/AnalysisTask.h"
@@ -30,6 +34,7 @@
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
 #include "PWGHF/Utils/utilsBfieldCCDB.h"
 #include "PWGHF/Utils/utilsTrkCandHf.h"
+#include "PWGHF/Utils/utilsMcGen.h"
 
 using namespace o2;
 using namespace o2::analysis;
@@ -448,20 +453,7 @@ struct HfCandidateCreatorB0Expressions {
       rowMcMatchRec(flag, origin, debug);
     } // rec
 
-    // Match generated particles.
-    for (const auto& particle : mcParticles) {
-      flag = 0;
-      origin = 0;
-      // B0 → D- π+
-      if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kB0, std::array{-static_cast<int>(Pdg::kDPlus), +kPiPlus}, true)) {
-        // D- → π- K+ π-
-        auto candDMC = mcParticles.rawIteratorAt(particle.daughtersIds().front());
-        if (RecoDecay::isMatchedMCGen(mcParticles, candDMC, -static_cast<int>(Pdg::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) {
-          flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi);
-        }
-      }
-      rowMcMatchGen(flag, origin);
-    } // gen
+    hf_mc_gen::fillMcMatchGenB0(mcParticles, rowMcMatchGen); // gen
   }   // processMc
   PROCESS_SWITCH(HfCandidateCreatorB0Expressions, processMc, "Process MC", false);
 }; // struct

PWGHF/TableProducer/candidateCreatorBplus.cxx

@@ -37,6 +37,7 @@
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
 #include "PWGHF/Utils/utilsBfieldCCDB.h"
 #include "PWGHF/Utils/utilsTrkCandHf.h"
+#include "PWGHF/Utils/utilsMcGen.h"
 
 using namespace o2;
 using namespace o2::analysis;
@@ -365,7 +366,6 @@ struct HfCandidateCreatorBplusExpressions {
     int8_t signB = 0, signD0 = 0;
     int8_t flag = 0;
     int8_t origin = 0;
-    int kD0pdg = Pdg::kD0;
 
     // Match reconstructed candidates.
     // Spawned table can be used directly
@@ -386,31 +386,7 @@ struct HfCandidateCreatorBplusExpressions {
       }
       rowMcMatchRec(flag, origin);
     }
-
-    // Match generated particles.
-    for (const auto& particle : mcParticles) {
-      flag = 0;
-      origin = 0;
-      signB = 0;
-      signD0 = 0;
-      int indexGenD0 = -1;
-
-      // B± → D0bar(D0) π± → (K± π∓) π±
-      std::vector<int> arrayDaughterB;
-      if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kBPlus, std::array{-kD0pdg, +kPiPlus}, true, &signB, 1, &arrayDaughterB)) {
-        // D0(bar) → π± K∓
-        for (auto iD : arrayDaughterB) {
-          auto candDaughterMC = mcParticles.rawIteratorAt(iD);
-          if (std::abs(candDaughterMC.pdgCode()) == kD0pdg) {
-            indexGenD0 = RecoDecay::isMatchedMCGen(mcParticles, candDaughterMC, Pdg::kD0, std::array{-kKPlus, +kPiPlus}, true, &signD0, 1);
-          }
-        }
-        if (indexGenD0 > -1) {
-          flag = signB * (1 << hf_cand_bplus::DecayType::BplusToD0Pi);
-        }
-      }
-      rowMcMatchGen(flag, origin);
-    } // B candidate
+    hf_mc_gen::fillMcMatchGenBplus(mcParticles, rowMcMatchGen); // gen
   }   // process
   PROCESS_SWITCH(HfCandidateCreatorBplusExpressions, processMc, "Process MC", false);
 }; // struct

PWGHF/TableProducer/candidateCreatorMcGen.cxx

@@ -0,0 +1,76 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file candidateCreatorMcGen.cxx
+/// \brief McGen only selection of heavy-flavour particles
+///
+/// \author Nima Zardoshti, nima.zardoshti@cern.ch, CERN
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <TPDGCode.h>
+
+#include "CommonConstants/PhysicsConstants.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/runDataProcessing.h"
+#include "Framework/RunningWorkflowInfo.h"
+
+#include "PWGHF/DataModel/CandidateReconstructionTables.h"
+#include "PWGHF/Utils/utilsMcGen.h"
+
+using namespace o2;
+using namespace o2::analysis;
+using namespace o2::framework;
+
+/// Reconstruction of heavy-flavour 2-prong decay candidates
+struct HfCandidateCreatorMcGen {
+
+  Produces<aod::HfCand2ProngMcGen> rowMcMatchGen2Prong;
+  Produces<aod::HfCand3ProngMcGen> rowMcMatchGen3Prong;
+  Produces<aod::HfCandBplusMcGen> rowMcMatchGenBplus;
+  Produces<aod::HfCandB0McGen> rowMcMatchGenB0;
+  Configurable<bool> fill2Prong{"fill2Prong", false, "fill table for 2 prong candidates"};
+  Configurable<bool> fill3Prong{"fill3Prong", false, "fill table for 3 prong candidates"};
+  Configurable<bool> fillBplus{"fillBplus", false, "fill table for for B+ candidates"};
+  Configurable<bool> fillB0{"fillB0", false, "fill table for B0 candidates"};
+  Configurable<bool> rejectBackground2Prong{"rejectBackground2Prong", false, "Reject particles from PbPb background for 2 prong candidates"};
+  Configurable<bool> rejectBackground3Prong{"rejectBackground3Prong", false, "Reject particles from PbPb background for 3 prong candidates"};
+  Configurable<bool> createDplus{"createDplus", false, "Create D+ in 3 prong"};
+  Configurable<bool> createDs{"createDs", false, "Create Ds in 3 prong"};
+  Configurable<bool> createLc{"createLc", false, "Create Lc in 3 prong"};
+  Configurable<bool> createXic{"createXic", false, "Create Xic in 3 prong"};
+
+  void process(aod::McCollision const&,
+               aod::McParticles const& mcParticles)
+  {
+    if (fill2Prong) {
+      hf_mc_gen::fillMcMatchGen2Prong(mcParticles, rowMcMatchGen2Prong, rejectBackground2Prong);
+    }
+    if (fill3Prong) {
+      hf_mc_gen::fillMcMatchGen3Prong(mcParticles, rowMcMatchGen3Prong, rejectBackground3Prong, createDplus, createDs, createLc, createXic);
+    }
+    if (fillBplus) {
+      hf_mc_gen::fillMcMatchGenBplus(mcParticles, rowMcMatchGenBplus);
+    }
+    if (fillB0) {
+      hf_mc_gen::fillMcMatchGenB0(mcParticles, rowMcMatchGenB0);
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<HfCandidateCreatorMcGen>(cfgc)};
+}

PWGHF/TableProducer/derivedDataCreatorBplusToD0Pi.cxx

@@ -419,6 +419,13 @@ struct HfDerivedDataCreatorBplusToD0Pi {
     rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass);
   }
   PROCESS_SWITCH(HfDerivedDataCreatorBplusToD0Pi, processMcMlAll, "Process MC with ML", false);
+
+  void processMcGenOnly(TypeMcCollisions const& mcCollisions,
+                        MatchedGenCandidatesMc const& mcParticles)
+  {
+    rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass);
+  }
+  PROCESS_SWITCH(HfDerivedDataCreatorBplusToD0Pi, processMcGenOnly, "Process MC gen. only", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

PWGHF/TableProducer/derivedDataCreatorD0ToKPi.cxx

@@ -516,6 +516,13 @@ struct HfDerivedDataCreatorD0ToKPi {
     rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass);
   }
   PROCESS_SWITCH(HfDerivedDataCreatorD0ToKPi, processMcWithKFParticleMlAll, "Process MC with KFParticle and ML", false);
+
+  void processMcGenOnly(TypeMcCollisions const& mcCollisions,
+                        MatchedGenCandidatesMc const& mcParticles)
+  {
+    rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass);
+  }
+  PROCESS_SWITCH(HfDerivedDataCreatorD0ToKPi, processMcGenOnly, "Process MC gen. only", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

PWGHF/TableProducer/derivedDataCreatorLcToPKPi.cxx

@@ -387,6 +387,13 @@ struct HfDerivedDataCreatorLcToPKPi {
     rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass);
   }
   PROCESS_SWITCH(HfDerivedDataCreatorLcToPKPi, processMcMlAll, "Process MC with ML", false);
+
+  void processMcGenOnly(TypeMcCollisions const& mcCollisions,
+                        MatchedGenCandidatesMc const& mcParticles)
+  {
+    rowsCommon.processMcParticles(mcCollisions, mcParticlesPerMcCollision, mcParticles, mass);
+  }
+  PROCESS_SWITCH(HfDerivedDataCreatorLcToPKPi, processMcGenOnly, "Process MC gen. only", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)