PWGJE: Add task for MC generator studies and write EMCreadout status …

…to table

- Add task that extracts generated pt spectrum of a specific particle
- Add isemcreadout boolean to emcalcollisionmatching table

PWGJE/DataModel/EMCALMatchedCollisions.h

@@ -26,12 +26,13 @@ namespace o2::aod
 {
 namespace emcalcollisionmatch
 {
-DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! collisionID used as index for matched collisions
-DECLARE_SOA_COLUMN(Ambiguous, ambiguous, bool); //! boolean stating whether the collision is ambiguous (in a BC with multiple collisions)
+DECLARE_SOA_INDEX_COLUMN(Collision, collision);       //! collisionID used as index for matched collisions
+DECLARE_SOA_COLUMN(Ambiguous, ambiguous, bool);       //! boolean stating whether the collision is ambiguous (in a BC with multiple collisions)
+DECLARE_SOA_COLUMN(IsEMCReadout, isemcreadout, bool); //! boolean stating whether the EMCal was readout in that collision (based on whether the BC contains at least one cell)
 } // namespace emcalcollisionmatch
 
-DECLARE_SOA_TABLE(EMCALMatchedCollisions, "AOD", "EMCALMCS",                                           //!
-                  o2::soa::Index<>, emcalcollisionmatch::CollisionId, emcalcollisionmatch::Ambiguous); //
+DECLARE_SOA_TABLE(EMCALMatchedCollisions, "AOD", "EMCALMCS",                                                                              //!
+                  o2::soa::Index<>, emcalcollisionmatch::CollisionId, emcalcollisionmatch::Ambiguous, emcalcollisionmatch::IsEMCReadout); //
 
 using EMCALMatchedCollision = EMCALMatchedCollisions::iterator;
 

PWGJE/TableProducer/emcalCorrectionTask.cxx

@@ -223,6 +223,7 @@ struct EmcalCorrectionTask {
     int nBCsProcessed = 0;
     int nCellsProcessed = 0;
     std::unordered_map<uint64_t, int> numberCollsInBC; // Number of collisions mapped to the global BC index of all BCs
+    std::unordered_map<uint64_t, int> numberCellsInBC; // Number of cells mapped to the global BC index of all BCs to check whether EMCal was readout
     for (auto bc : bcs) {
       LOG(debug) << "Next BC";
       // Convert aod::Calo to o2::emcal::Cell which can be used with the clusterizer.
@@ -233,6 +234,7 @@ struct EmcalCorrectionTask {
       auto cellsInBC = cells.sliceBy(cellsPerFoundBC, bc.globalIndex());
 
       numberCollsInBC.insert(std::pair<uint64_t, int>(bc.globalIndex(), collisionsInFoundBC.size()));
+      numberCellsInBC.insert(std::pair<uint64_t, int>(bc.globalIndex(), cellsInBC.size()));
 
       if (!cellsInBC.size()) {
         LOG(debug) << "No cells found for BC";
@@ -316,9 +318,10 @@ struct EmcalCorrectionTask {
     // Loop through all collisions and fill emcalcollisionmatch with a boolean stating, whether the collision was ambiguous (not the only collision in its BC)
     for (const auto& collision : collisions) {
       auto globalbcid = collision.foundBC_as<bcEvSels>().globalIndex();
-      auto found = numberCollsInBC.find(globalbcid);
-      if (found != numberCollsInBC.end()) {
-        emcalcollisionmatch(collision.globalIndex(), found->second != 1);
+      auto foundColls = numberCollsInBC.find(globalbcid);
+      auto foundCells = numberCellsInBC.find(globalbcid);
+      if (foundColls != numberCollsInBC.end() && foundCells != numberCellsInBC.end()) {
+        emcalcollisionmatch(collision.globalIndex(), foundColls->second != 1, foundCells->second > 0);
       } else {
         LOG(warning) << "BC not found in map of number of collisions.";
       }

PWGJE/Tasks/CMakeLists.txt

@@ -34,6 +34,10 @@ o2physics_add_dpl_workflow(emc-tmmonitor
                     SOURCES emctmmonitor.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
+o2physics_add_dpl_workflow(mc-generator-studies
+                    SOURCES mcgeneratorstudies.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
 
 if(FastJet_FOUND)
     o2physics_add_dpl_workflow(jet-substructure

PWGJE/Tasks/mcgeneratorstudies.cxx

@@ -0,0 +1,88 @@
+// Copyright 2019-2024 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.
+
+// Task that produces the generated pT spectrum of a given particle for MC studies
+//
+/// \author Nicolas Strangmann <[email protected]>, Goethe University Frankfurt / Oak Ridge National Laoratory
+
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/ASoA.h"
+#include "Framework/HistogramRegistry.h"
+
+#include "PWGJE/DataModel/EMCALMatchedCollisions.h"
+
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/DataModel/EventSelection.h"
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+using MyMCCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::EMCALMatchedCollisions>;
+
+struct MCGeneratorStudies {
+  HistogramRegistry mHistManager{"MCGeneratorStudyHistograms"};
+
+  Configurable<float> mVertexCut{"vertexCut", 10.f, "apply z-vertex cut with value in cm"};
+  Configurable<double> mSelectedParticleCode{"particlePDGCode", 111, "PDG code of the particle to be investigated"};
+
+  Filter collisionFilter = (aod::collision::posZ > -mVertexCut) && (aod::collision::posZ < mVertexCut);
+  Filter mcParticleFilter = aod::mcparticle::pdgCode == mSelectedParticleCode;
+
+  void init(InitContext const&)
+  {
+    AxisSpec pTAxis{250, 0., 25., "#it{p}_{T} (GeV/#it{c})"};
+
+    mHistManager.add("hCollisionCounter", "Number of collisions after event cuts", HistType::kTH1F, {{5, 0.5, 5.5}});
+    mHistManager.add("hpT_all", "All collisions", HistType::kTH1F, {pTAxis});
+    mHistManager.add("hpT_T0Triggered", "T0 triggered collisions", HistType::kTH1F, {pTAxis});
+    mHistManager.add("hpT_T0Triggered_EMCReadout", "T0 triggered and EMC readout collisions", HistType::kTH1F, {pTAxis});
+    mHistManager.add("hpT_T0Triggered_EMCReadout_Unique", "Unique T0 triggered and EMC readout collisions", HistType::kTH1F, {pTAxis});
+  }
+
+  void process(soa::Filtered<MyMCCollisions>::iterator const& collision, soa::Filtered<aod::McParticles> const& mcParticles)
+  {
+    bool isT0Triggered = collision.sel8();
+    bool isEMCReadout = collision.isemcreadout();
+    bool isUniqueCollision = !collision.ambiguous();
+
+    mHistManager.fill(HIST("hCollisionCounter"), 1);
+    if (isT0Triggered) {
+      mHistManager.fill(HIST("hCollisionCounter"), 2);
+      // mHistManager.fill<aod::mcparticle::Pt>(HIST("hpT_T0Triggered"), mcParticles, aod::evsel::sel8 == true);
+      if (isEMCReadout) {
+        mHistManager.fill(HIST("hCollisionCounter"), 3);
+        if (isUniqueCollision) {
+          mHistManager.fill(HIST("hCollisionCounter"), 4);
+        }
+      }
+    }
+    for (auto& mcParticle : mcParticles) {
+      mHistManager.fill(HIST("hpT_all"), mcParticle.pt());
+      if (isT0Triggered) {
+        mHistManager.fill(HIST("hpT_T0Triggered"), mcParticle.pt());
+        if (isEMCReadout) {
+          mHistManager.fill(HIST("hpT_T0Triggered_EMCReadout"), mcParticle.pt());
+          if (isUniqueCollision) {
+            mHistManager.fill(HIST("hpT_T0Triggered_EMCReadout_Unique"), mcParticle.pt());
+          }
+        }
+      }
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<MCGeneratorStudies>(cfgc, TaskName{"mc-generator-studies"})};
+}