diff --git a/PWGEM/PhotonMeson/DataModel/gammaTables.h b/PWGEM/PhotonMeson/DataModel/gammaTables.h index 0851203654e..e1ca0cc8d54 100644 --- a/PWGEM/PhotonMeson/DataModel/gammaTables.h +++ b/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_ diff --git a/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt b/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt index dc400c3ba95..1e9c99c1cd5 100644 --- a/PWGEM/PhotonMeson/TableProducer/CMakeLists.txt +++ b/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) diff --git a/PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx b/PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx new file mode 100644 index 00000000000..0561e60fc22 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +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 tableNonLinV0s; + Produces tableNonLinClusters; + + Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3)"}; + + HistogramRegistry historeg{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + using EMCalPhotons = soa::Join; + using PcmPhotons = soa::Join; + + using Colls = soa::Join; + + 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 + 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 + 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 + 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(cfgc)}; + return workflow; +} diff --git a/PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h b/PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h index 0fc655d8858..846168d774c 100644 --- a/PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h +++ b/PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h @@ -39,12 +39,23 @@ const std::vector 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 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}});