[PWGEM] PhotonMeson: Add NonLin tables and producer

PWGEM/PhotonMeson/DataModel/gammaTables.h

@@ -682,6 +682,20 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Cell, cell, int, Calos, "");                 //! r
 DECLARE_SOA_TABLE(SkimEMCCells, "AOD", "SKIMEMCCELLS",                        //! table of link between skimmed EMCal clusters and their cells
                   o2::soa::Index<>, caloextra::ClusterId, caloextra::CellId); //!
 using SkimEMCCell = SkimEMCCells::iterator;
+
+namespace nonlin
+{
+DECLARE_SOA_COLUMN(CorrE, corrE, float);   //! non lin corrected energy
+DECLARE_SOA_COLUMN(CorrP, corrP, float);   //! non lin corrected momentum
+DECLARE_SOA_COLUMN(CorrPt, corrPt, float); //! non lin corrected transverse momentum
+} // namespace nonlin
+
+DECLARE_SOA_TABLE(NonLinV0s, "AOD", "NONLINV0", //! table of non lin corrected values for V0 Photons (so far only pt)
+                  nonlin::CorrPt);              //!
+
+DECLARE_SOA_TABLE(NonLinEmcClusters, "AOD", "NONLINEMCCLUSTER", //! table of non lin corrected values for EMCal Photons (so far only E and pT)
+                  nonlin::CorrE, nonlin::CorrPt);               //!
+
 } // namespace o2::aod
 
 #endif // PWGEM_PHOTONMESON_DATAMODEL_GAMMATABLES_H_

PWGEM/PhotonMeson/TableProducer/CMakeLists.txt

@@ -63,3 +63,8 @@ o2physics_add_dpl_workflow(skimmer-dalitz-ee
                     SOURCES skimmerDalitzEE.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(non-lin-producer
+                    SOURCES nonLinProducer.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore
+                    COMPONENT_NAME Analysis)

PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx

@@ -0,0 +1,206 @@
+// 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 nonLinProducer.cxx
+/// \brief produces nonLin tables for PCM and EMCal
+/// \author marvin.hemmer@cern.ch
+/// dependencies: skimmer-gamma-calo
+
+#include "PWGEM/PhotonMeson/Core/EMNonLin.h"
+#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h"
+#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
+#include "PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h"
+
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/Configurable.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/InitContext.h>
+#include <Framework/OutputObjHeader.h>
+#include <Framework/runDataProcessing.h>
+
+#include <cstdint>
+#include <vector>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::aod::emcdownscaling;
+using namespace o2::pwgem::nonlin;
+
+struct NonLinProducer {
+
+  enum CentralityEstimator : uint8_t {
+    None = 0,
+    CFT0A = 1,
+    CFT0C,
+    CFT0M,
+    NCentralityEstimators
+  };
+
+  Produces<aod::NonLinV0s> tableNonLinV0s;
+  Produces<aod::NonLinEmcClusters> tableNonLinClusters;
+
+  Configurable<int> centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3)"};
+
+  HistogramRegistry historeg{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
+
+  using EMCalPhotons = soa::Join<aod::EMCEMEventIds, aod::MinClusters>;
+  using PcmPhotons = soa::Join<aod::V0PhotonsKF, aod::V0KFEMEventIds>;
+
+  using Colls = soa::Join<aod::EMEvents, aod::EMEventsCent>;
+
+  EMNonLin emNonLinEMC;
+  EMNonLin emNonLinPCM;
+
+  void init(o2::framework::InitContext&)
+  {
+    historeg.add("QA/EMC/EIn", "Energy of clusters before cuts", gHistoSpecClusterE);
+    historeg.add("QA/EMC/EOut", "Energy of clusters after cuts", gHistoSpecClusterE);
+
+    historeg.add("QA/EMC/PtIn", "Energy of clusters before cuts", gHistoSpecPt);
+    historeg.add("QA/EMC/PtOut", "Energy of clusters after cuts", gHistoSpecPt);
+
+    historeg.add("QA/PCM/PtIn", "Energy of clusters before cuts", gHistoSpecPt);
+    historeg.add("QA/PCM/PtOut", "Energy of clusters after cuts", gHistoSpecPt);
+  }
+
+  /// Get the centrality
+  /// \param collision is the collision with the centrality information
+  template <o2::soa::is_iterator TCollision>
+  float getCentrality(TCollision const& collision)
+  {
+    float cent = -999.;
+    switch (centEstimator) {
+      case CentralityEstimator::CFT0M:
+        cent = collision.centFT0M();
+        break;
+      case CentralityEstimator::CFT0A:
+        cent = collision.centFT0A();
+        break;
+      case CentralityEstimator::CFT0C:
+        cent = collision.centFT0C();
+        break;
+      default:
+        LOG(warning) << "Centrality estimator not valid. Possible values are T0M, T0A, T0C. Fallback to T0C";
+        cent = collision.centFT0C();
+        break;
+    }
+    return cent;
+  }
+
+  template <o2::soa::is_table TClusters, o2::soa::is_iterator TCollisio>
+  void runEMC(TClusters const& clusters, TCollisio& collision)
+  {
+    float nonLinE = 0.f;
+    float nonLinPt = 0.f;
+
+    float nonLinFactor = 1.f;
+
+    int32_t collIndex = collision.globalIndex();
+    for (const auto& cluster : clusters) {
+
+      // check that we are at the correct collision
+      if (cluster.emeventId() != collIndex) {
+        collIndex = cluster.emeventId();
+        collision.setCursor(collIndex);
+      }
+
+      // fill before non lin histograms
+      historeg.fill(HIST("QA/EMC/EIn"), cluster.e());
+      historeg.fill(HIST("QA/EMC/PtIn"), cluster.pt());
+
+      // get NonLin factor from class dependent on the centrality
+      nonLinFactor = emNonLinEMC.getCorrectionFactor(cluster.e(), o2::pwgem::nonlin::EMNonLin::PhotonType::kEMC, getCentrality(collision));
+
+      nonLinE = nonLinFactor * cluster.e();
+      nonLinPt = nonLinFactor * cluster.pt();
+
+      // fill after non lin histograms
+      historeg.fill(HIST("QA/EMC/EIn"), nonLinE);
+      historeg.fill(HIST("QA/EMC/PtIn"), nonLinPt);
+
+      tableNonLinClusters(nonLinE, nonLinPt);
+    }
+  }
+
+  template <o2::soa::is_table TV0, o2::soa::is_iterator TCollisio>
+  void runPCM(TV0 const& v0s, TCollisio& collision)
+  {
+    float nonLinPt = 0.f;
+
+    float nonLinFactor = 1.f;
+
+    int32_t collIndex = collision.globalIndex();
+    for (const auto& v0 : v0s) {
+
+      // check that we are at the correct collision
+      if (v0.emeventId() != collIndex) {
+        collIndex = v0.emeventId();
+        collision.setCursor(collIndex);
+      }
+
+      // fill before non lin histograms
+      historeg.fill(HIST("QA/PCM/PtIn"), v0.pt());
+
+      // get NonLin factor from class dependent on the centrality
+      nonLinFactor = emNonLinEMC.getCorrectionFactor(v0.pt(), o2::pwgem::nonlin::EMNonLin::PhotonType::kPCM, getCentrality(collision));
+
+      nonLinPt = nonLinFactor * v0.pt();
+
+      // fill after non lin histograms
+      historeg.fill(HIST("QA/PCM/PtIn"), nonLinPt);
+
+      tableNonLinV0s(nonLinPt);
+    }
+  }
+
+  void processEMC(Colls const& collisions, EMCalPhotons const& emcclusters)
+  {
+
+    if (emcclusters.size() == 0) {
+      return;
+    }
+
+    auto collision = collisions.begin();
+    runEMC(emcclusters, collision);
+  }
+  PROCESS_SWITCH(NonLinProducer, processEMC, "Create Non Lin table for EMC.", false);
+
+  void processPCM(Colls const& collisions, PcmPhotons const& pcmPhotons)
+  {
+    if (pcmPhotons.size() == 0) {
+      return;
+    }
+    auto collision = collisions.begin();
+    runPCM(pcmPhotons, collision);
+  }
+  PROCESS_SWITCH(NonLinProducer, processPCM, "Create Non Lin table for PCM.", false);
+
+  void processEMCDummy(EMCalPhotons::iterator const& emccluster)
+  {
+    tableNonLinClusters(emccluster.e(), emccluster.pt());
+  }
+  PROCESS_SWITCH(NonLinProducer, processEMCDummy, "Create dummy Non Lin table for EMC.", true);
+
+  void processPCMDummy(PcmPhotons::iterator const& pcmPhoton)
+  {
+    tableNonLinV0s(pcmPhoton.pt());
+  }
+  PROCESS_SWITCH(NonLinProducer, processPCMDummy, "Createdumy  Non Lin table for PCM.", true);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  WorkflowSpec workflow{adaptAnalysisTask<NonLinProducer>(cfgc)};
+  return workflow;
+}

PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h

@@ -39,12 +39,23 @@ const std::vector<double> eClusBins = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8
                                        100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200};
 o2::framework::AxisSpec const gAxisEClus = {eClusBins, "#it{E}_{clus} (GeV)"};
 
+const std::vector<double> pTBins = {0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,
+                                    0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
+                                    1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8,
+                                    2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8,
+                                    3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8,
+                                    4.9, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9,
+                                    9.5, 10, 11, 12, 13, 14, 15, 16, 17, 18,
+                                    19, 20};
+o2::framework::AxisSpec const gAxisPt = {pTBins, "#it{p}_{T} (GeV/#it{c})"};
+
 o2::framework::HistogramConfigSpec const gHistoSpecClusterTMdEtadPhi({o2::framework::HistType::kTH2F, {gAxisdEta, gAxisdPhi}});
 o2::framework::HistogramConfigSpec const gHistoSpecClusterTMdEtaMT({o2::framework::HistType::kTH2F, {gAxisdEta, {10, 0.5, 10.5}}});
 o2::framework::HistogramConfigSpec const gHistoSpecClusterTMdPhiMT({o2::framework::HistType::kTH2F, {gAxisdPhi, {10, 0.5, 10.5}}});
 o2::framework::HistogramConfigSpec const gHistoSpecClusterTMdRMT({o2::framework::HistType::kTH2F, {gAxisdR, {10, 0.5, 10.5}}});
 o2::framework::HistogramConfigSpec const gHistoSpecEtaPhi({o2::framework::HistType::kTH2F, {gAxisEta, gAxisPhi}});
 o2::framework::HistogramConfigSpec const gHistoSpecClusterE({o2::framework::HistType::kTH1F, {gAxisEClus}});
+o2::framework::HistogramConfigSpec const gHistoSpecPt({o2::framework::HistType::kTH1F, {gAxisPt}});
 o2::framework::HistogramConfigSpec const gHistoSpecClusterECuts({o2::framework::HistType::kTH2F, {gAxisEClus, {64, -0.5, 63.5}}});
 o2::framework::HistogramConfigSpec const gHistoSpecTMEoverP({o2::framework::HistType::kTH2D, {gAxisEClus, gAxisEoverP}});