[PWGEM/Dilepton] add status code in mc particle table

PWGEM/Dilepton/DataModel/dileptonTables.h

@@ -64,15 +64,15 @@
 namespace emevent
 {
 DECLARE_SOA_COLUMN(CollisionId, collisionId, int);
 DECLARE_SOA_BITMAP_COLUMN(SWTAliasTmp, swtaliastmp, 16);             //! Bitmask of fired trigger aliases (see above for definitions) to be join to aod::Collisions for skimming
 DECLARE_SOA_BITMAP_COLUMN(SWTAlias, swtalias, 16);                   //! Bitmask of fired trigger aliases (see above for definitions) to be join to aod::EMEvents for analysis
 DECLARE_SOA_COLUMN(NInspectedTVX, nInspectedTVX, uint64_t);          //! the number of inspected TVX bcs per run
 DECLARE_SOA_COLUMN(NScalars, nScalers, std::vector<uint64_t>);       //! the number of triggered bcs before down scaling per run
 DECLARE_SOA_COLUMN(NSelections, nSelections, std::vector<uint64_t>); //! the number of triggered bcs after down scaling per run
 DECLARE_SOA_BITMAP_COLUMN(IsAnalyzed, isAnalyzed, 16);
 DECLARE_SOA_BITMAP_COLUMN(IsAnalyzedToI, isAnalyzedToI, 16);
 DECLARE_SOA_COLUMN(NeeULS, neeuls, int);
 DECLARE_SOA_COLUMN(NeeLSpp, neelspp, int);
 DECLARE_SOA_COLUMN(NeeLSmm, neelsmm, int);
 DECLARE_SOA_COLUMN(Bz, bz, float);                                          //! kG
 DECLARE_SOA_COLUMN(Q2xFT0M, q2xft0m, float);                                //! Qx for 2nd harmonics in FT0M
@@ -301,7 +301,7 @@
 } // namespace emmcparticle
 
 // This table contains all MC truth tracks
-DECLARE_SOA_TABLE_FULL(EMMCParticles, "EMMCParticles", "AOD", "EMMCPARTICLE", //!  MC track information (on disk)
+DECLARE_SOA_TABLE_FULL(EMMCParticles_000, "EMMCParticles", "AOD", "EMMCPARTICLE", //!  MC track information (on disk)
                        o2::soa::Index<>, emmcparticle::EMMCEventId,
                        mcparticle::PdgCode, mcparticle::Flags,
                        emmcparticle::MothersIds, emmcparticle::DaughtersIds,
@@ -312,13 +312,33 @@
                        emmcparticle::Pt<mcparticle::Px, mcparticle::Py>,
                        emmcparticle::Eta<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
                        emmcparticle::Phi<mcparticle::Px, mcparticle::Py>,
-
                        emmcparticle::P<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
                        emmcparticle::Y<mcparticle::Pz, mcparticle::E>,
                        mcparticle::ProducedByGenerator<mcparticle::Flags>,
                        mcparticle::FromBackgroundEvent<mcparticle::Flags>,
                        mcparticle::IsPhysicalPrimary<mcparticle::Flags>);
 
+DECLARE_SOA_TABLE_VERSIONED(EMMCParticles_001, "AOD", "EMMCPARTICLE", 1, //!  MC track information (on disk)
+                            o2::soa::Index<>, emmcparticle::EMMCEventId,
+                            mcparticle::PdgCode, mcparticle::Flags, mcparticle::StatusCode,
+                            emmcparticle::MothersIds, emmcparticle::DaughtersIds,
+                            mcparticle::Px, mcparticle::Py, mcparticle::Pz, mcparticle::E,
+                            mcparticle::Vx, mcparticle::Vy, mcparticle::Vz,
+
+                            // dynamic column
+                            emmcparticle::Pt<mcparticle::Px, mcparticle::Py>,
+                            emmcparticle::Eta<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
+                            emmcparticle::Phi<mcparticle::Px, mcparticle::Py>,
+                            emmcparticle::P<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
+                            emmcparticle::Y<mcparticle::Pz, mcparticle::E>,
+                            mcparticle::ProducedByGenerator<mcparticle::Flags>,
+                            mcparticle::FromBackgroundEvent<mcparticle::Flags>,
+                            mcparticle::IsPhysicalPrimary<mcparticle::Flags>,
+                            mcparticle::GetGenStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
+                            mcparticle::GetHepMCStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
+                            mcparticle::GetProcess<mcparticle::Flags, mcparticle::StatusCode>);
+
+using EMMCParticles = EMMCParticles_001;
 using EMMCParticle = EMMCParticles::iterator;
 
 namespace emmcgenvectormeson
@@ -437,7 +457,7 @@
 DECLARE_SOA_DYNAMIC_COLUMN(Tgl, tgl, [](float eta) -> float { return std::tan(M_PI_2 - 2 * std::atan(std::exp(-eta))); });
 DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITS, meanClusterSizeITS, [](uint32_t itsClusterSizes) -> float {
   int total_cluster_size = 0, nl = 0;
   for (unsigned int layer = 0; layer < 7; layer++) {
     int cluster_size_per_layer = (itsClusterSizes >> (layer * 4)) & 0xf;
     if (cluster_size_per_layer > 0) {
       nl++;
@@ -452,7 +472,7 @@
 });
 DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSib, meanClusterSizeITSib, [](uint32_t itsClusterSizes) -> float {
   int total_cluster_size = 0, nl = 0;
   for (unsigned int layer = 0; layer < 3; layer++) {
     int cluster_size_per_layer = (itsClusterSizes >> (layer * 4)) & 0xf;
     if (cluster_size_per_layer > 0) {
       nl++;
@@ -467,7 +487,7 @@
 });
 DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSob, meanClusterSizeITSob, [](uint32_t itsClusterSizes) -> float {
   int total_cluster_size = 0, nl = 0;
   for (unsigned int layer = 3; layer < 7; layer++) {
     int cluster_size_per_layer = (itsClusterSizes >> (layer * 4)) & 0xf;
     if (cluster_size_per_layer > 0) {
       nl++;
@@ -735,7 +755,7 @@
 DECLARE_SOA_DYNAMIC_COLUMN(NClustersMFT, nClustersMFT, //! Number of MFT clusters
                            [](uint64_t mftClusterSizesAndTrackFlags) -> uint8_t {
                              uint8_t nClusters = 0;
                              for (int layer = 0; layer < 10; layer++) {
                                if ((mftClusterSizesAndTrackFlags >> (layer * 6)) & 0x3F) {
                                  nClusters++;
                                }
@@ -745,7 +765,7 @@
 DECLARE_SOA_DYNAMIC_COLUMN(MFTClusterMap, mftClusterMap, //! MFT cluster map, one bit per a layer, starting from the innermost
                            [](uint64_t mftClusterSizesAndTrackFlags) -> uint16_t {
                              uint16_t clmap = 0;
                              for (unsigned int layer = 0; layer < 10; layer++) {
                                if ((mftClusterSizesAndTrackFlags >> (layer * 6)) & 0x3f) {
                                  clmap |= (1 << layer);
                                }

PWGEM/Dilepton/TableProducer/associateMCinfoDilepton.cxx

@@ -625,7 +625,7 @@ struct AssociateMCInfoDilepton {
         }
       }
 
-      emmcparticles(fEventIdx.find(oldLabel)->second, mctrack.pdgCode(), mctrack.flags(),
+      emmcparticles(fEventIdx.find(oldLabel)->second, mctrack.pdgCode(), mctrack.flags(), mctrack.statusCode(),
                     mothers, daughters,
                     mctrack.px(), mctrack.py(), mctrack.pz(), mctrack.e(),
                     mctrack.vx(), mctrack.vy(), mctrack.vz());

PWGEM/Dilepton/Tasks/Converters/CMakeLists.txt

@@ -50,3 +50,8 @@ o2physics_add_dpl_workflow(muon-converter1
                   PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
                   COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(mcparticle-converter1
+                  SOURCES mcParticleConverter1.cxx
+                  PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
+                  COMPONENT_NAME Analysis)
+

PWGEM/Dilepton/Tasks/Converters/mcParticleConverter1.cxx

@@ -0,0 +1,61 @@
+// 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.
+//
+// ========================
+//
+// This code produces emmctable table 001 from 000.
+//    Please write to: daiki.sekihata@cern.ch
+
+#include "PWGEM/Dilepton/DataModel/dileptonTables.h"
+
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+
+#include <vector>
+
+using namespace o2;
+using namespace o2::aod;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::soa;
+
+struct mcParticleConverter1 {
+  Produces<aod::EMMCParticles_001> mcParticle_001;
+
+  void process(aod::EMMCParticles_000 const& mcParticles)
+  {
+    for (const auto& mcParticle : mcParticles) {
+      // LOGF(info, "mcParticles.emmceventId() = %d, mcParticle.mothersIds().size() = %d, mcParticle.daughtersIds().size() = %d", mcParticle.emmceventId(), mcParticle.mothersIds().size(), mcParticle.daughtersIds().size());
+
+      std::vector<int> mothersIds;
+      for (const auto& id : mcParticle.mothersIds()) {
+        mothersIds.emplace_back(id);
+      }
+
+      std::vector<int> daughtersIds;
+      for (const auto& id : mcParticle.daughtersIds()) {
+        daughtersIds.emplace_back(id);
+      }
+
+      mcParticle_001(
+        mcParticle.emmceventId(), mcParticle.pdgCode(), mcParticle.flags(), 0,
+        mothersIds, daughtersIds,
+        mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e(),
+        mcParticle.vx(), mcParticle.vy(), mcParticle.vz());
+    } // end of mc particle loop
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<mcParticleConverter1>(cfgc, TaskName{"mcparticle-converter1"})};
+}

PWGEM/Dilepton/Utils/MCUtilities.h

@@ -317,7 +317,7 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp
       auto mp = mcparticles.iteratorAt(motherid1);
       mothers_id1.emplace_back(motherid1);
       mothers_pdg1.emplace_back(mp.pdgCode());
-
+      // LOGF(info, "mp1.globalIndex() = %d, mp1.pdgCode() = %d, mp1.getGenStatusCode() = %d", mp.globalIndex(), mp.pdgCode(), mp.getGenStatusCode());
       if (mp.has_mothers()) {
         motherid1 = mp.mothersIds()[0];
       } else {
@@ -337,6 +337,7 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp
       auto mp = mcparticles.iteratorAt(motherid2);
       mothers_id2.emplace_back(motherid2);
       mothers_pdg2.emplace_back(mp.pdgCode());
+      // LOGF(info, "mp2.globalIndex() = %d, mp2.pdgCode() = %d, mp2.getGenStatusCode() = %d", mp.globalIndex(), mp.pdgCode(), mp.getGenStatusCode());
 
       if (mp.has_mothers()) {
         motherid2 = mp.mothersIds()[0];
@@ -348,6 +349,8 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp
     }
   }
 
+  // require correlation between q-qbar. (not q-q) // need statusCode
+
   bool is_direct_from_b1 = IsFromBeauty(p1, mcparticles) > 0 && IsFromCharm(p1, mcparticles) < 0;
   bool is_direct_from_b2 = IsFromBeauty(p2, mcparticles) > 0 && IsFromCharm(p2, mcparticles) < 0;
   bool is_prompt_c1 = IsFromBeauty(p1, mcparticles) < 0 && IsFromCharm(p1, mcparticles) > 0;

PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

@@ -16,16 +16,17 @@
 /// \author Daiki Sekihata (daiki.sekihata@cern.ch)
 ///
 
-#include <map>
-#include <vector>
+#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
+#include "PWGEM/PhotonMeson/Utils/MCUtilities.h"
 
-#include "Framework/runDataProcessing.h"
-#include "Framework/AnalysisTask.h"
-#include "Framework/AnalysisDataModel.h"
 #include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
 #include "ReconstructionDataFormats/Track.h"
-#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
-#include "PWGEM/PhotonMeson/Utils/MCUtilities.h"
+
+#include <map>
+#include <vector>
 
 using namespace o2;
 using namespace o2::framework;
@@ -434,7 +435,7 @@ struct AssociateMCInfoPhoton {
         }
       }
 
-      emmcparticles(fEventIdx.find(oldLabel)->second, mcParticle.pdgCode(), mcParticle.flags(),
+      emmcparticles(fEventIdx.find(oldLabel)->second, mcParticle.pdgCode(), mcParticle.flags(), mcParticle.statusCode(),
                     mothers, daughters,
                     mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e(),
                     mcParticle.vx(), mcParticle.vy(), mcParticle.vz());