Skip to content

Commit a395f40

Browse files
romainschotterromainschotter
authored
[PWGDQ] Add quarkonia to hyperon-antihyperon analysis tasks (#8342)
1 parent 8517b9a commit a395f40

File tree

5 files changed

+2280
-0
lines changed

5 files changed

+2280
-0
lines changed

‎PWGDQ/DataModel/ReducedInfoTables.h‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1069,6 +1069,30 @@ DECLARE_SOA_TABLE(ReducedMuonsDca, "AOD", "RTMUONDCA",
10691069
muondca::Pz);
10701070

10711071
using ReducedMuonDca = ReducedMuonsDca::iterator;
1072+
1073+
//______________________________________________________
1074+
namespace generatedquarkoniamc
1075+
{
1076+
//______________________________________________________
1077+
// Binned content for generated particles: derived data
1078+
DECLARE_SOA_COLUMN(GeneratedEtaC1S, generatedEtaC1S, std::vector<uint32_t>); //! Eta(1S) binned generated data
1079+
DECLARE_SOA_COLUMN(GeneratedJPsi, generatedJPsi, std::vector<uint32_t>); //! J/Psi binned generated data
1080+
DECLARE_SOA_COLUMN(GeneratedChiC0, generatedChiC0, std::vector<uint32_t>); //! ChiC0(1P) binned generated data
1081+
DECLARE_SOA_COLUMN(GeneratedChiC1, generatedChiC1, std::vector<uint32_t>); //! ChiC1(1P) binned generated data
1082+
DECLARE_SOA_COLUMN(GeneratedHC, generatedHC, std::vector<uint32_t>); //! hC binned generated data
1083+
DECLARE_SOA_COLUMN(GeneratedChiC2, generatedChiC2, std::vector<uint32_t>); //! ChiC2(1P) binned generated data
1084+
DECLARE_SOA_COLUMN(GeneratedEtaC2S, generatedEtaC2S, std::vector<uint32_t>); //! EtaC(2S) binned generated data
1085+
DECLARE_SOA_COLUMN(GeneratedPsi2S, generatedPsi2S, std::vector<uint32_t>); //! Psi(2S) binned generated data
1086+
} // namespace generatedquarkoniamc
1087+
1088+
DECLARE_SOA_TABLE(GeEtaC1S, "AOD", "GEETAC1S", generatedquarkoniamc::GeneratedEtaC1S);
1089+
DECLARE_SOA_TABLE(GeJPsi, "AOD", "GEJPSI", generatedquarkoniamc::GeneratedJPsi);
1090+
DECLARE_SOA_TABLE(GeChiC0, "AOD", "GECHIC0", generatedquarkoniamc::GeneratedChiC0);
1091+
DECLARE_SOA_TABLE(GeChiC1, "AOD", "GECHIC1", generatedquarkoniamc::GeneratedChiC1);
1092+
DECLARE_SOA_TABLE(GeHC, "AOD", "GEHC", generatedquarkoniamc::GeneratedHC);
1093+
DECLARE_SOA_TABLE(GeChiC2, "AOD", "GECHIC2", generatedquarkoniamc::GeneratedChiC2);
1094+
DECLARE_SOA_TABLE(GeEtaC2S, "AOD", "GEETAC2S", generatedquarkoniamc::GeneratedEtaC2S);
1095+
DECLARE_SOA_TABLE(GePsi2S, "AOD", "GEPSI2S", generatedquarkoniamc::GeneratedPsi2S);
10721096
} // namespace o2::aod
10731097

10741098
#endif // PWGDQ_DATAMODEL_REDUCEDINFOTABLES_H_

‎PWGDQ/TableProducer/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,3 +38,8 @@ o2physics_add_dpl_workflow(table-maker-muon-mch-trk-eff
3838
SOURCES tableMakerMuonMchTrkEfficiency.cxx
3939
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore
4040
COMPONENT_NAME Analysis)
41+
42+
o2physics_add_dpl_workflow(generated-quarkonia-mc
43+
SOURCES generatedQuarkoniaMC.cxx
44+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase O2Physics::PWGDQCore
45+
COMPONENT_NAME Analysis)

‎PWGDQ/TableProducer/generatedQuarkoniaMC.cxx‎

Lines changed: 242 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,242 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
//
12+
//__________________________________________________
13+
// this task provides produces histograms containing
14+
// the number of generated quarkonia per unit of
15+
// percentile and per unit of pT
16+
// It is meant to help with providing auxiliary information
17+
// when dealing with derived data.
18+
19+
#include <cmath>
20+
#include <array>
21+
#include <cstdlib>
22+
#include <map>
23+
#include <iterator>
24+
#include <utility>
25+
26+
#include "Framework/runDataProcessing.h"
27+
#include "Framework/RunningWorkflowInfo.h"
28+
#include "Framework/AnalysisTask.h"
29+
#include "Framework/AnalysisDataModel.h"
30+
#include "Framework/ASoAHelpers.h"
31+
#include "DCAFitter/DCAFitterN.h"
32+
#include "ReconstructionDataFormats/Track.h"
33+
#include "Common/Core/RecoDecay.h"
34+
#include "Common/Core/trackUtilities.h"
35+
#include "PWGLF/DataModel/LFStrangenessTables.h"
36+
#include "PWGLF/DataModel/LFStrangenessPIDTables.h"
37+
#include "PWGLF/DataModel/LFParticleIdentification.h"
38+
#include "PWGDQ/DataModel/ReducedInfoTables.h"
39+
#include "Common/Core/TrackSelection.h"
40+
#include "Common/DataModel/TrackSelectionTables.h"
41+
#include "DetectorsBase/Propagator.h"
42+
#include "DetectorsBase/GeometryManager.h"
43+
#include "DataFormatsParameters/GRPObject.h"
44+
#include "DataFormatsParameters/GRPMagField.h"
45+
#include "CCDB/BasicCCDBManager.h"
46+
#include "CommonConstants/PhysicsConstants.h"
47+
#include "Common/TableProducer/PID/pidTOFBase.h"
48+
#include "Common/DataModel/PIDResponse.h"
49+
#include "Common/DataModel/Qvectors.h"
50+
#include "Framework/StaticFor.h"
51+
#include "Common/DataModel/McCollisionExtra.h"
52+
#include "PWGLF/DataModel/EPCalibrationTables.h"
53+
54+
using namespace o2;
55+
using namespace o2::framework;
56+
using namespace o2::framework::expressions;
57+
using std::array;
58+
59+
// simple bit checkers
60+
#define bitset(var, nbit) ((var) |= (1 << (nbit)))
61+
#define bitcheck(var, nbit) ((var) & (1 << (nbit)))
62+
63+
struct generatedQuarkoniaMC {
64+
SliceCache cache;
65+
//__________________________________________________
66+
// Generated binned data
67+
// this is a hack while the system does not do better
68+
Produces<aod::GeEtaC1S> geEtaC1S;
69+
Produces<aod::GeJPsi> geJPsi;
70+
Produces<aod::GeChiC0> geChiC0;
71+
Produces<aod::GeChiC1> geChiC1;
72+
Produces<aod::GeHC> geHC;
73+
Produces<aod::GeChiC2> geChiC2;
74+
Produces<aod::GeEtaC2S> geEtaC2S;
75+
Produces<aod::GePsi2S> gePsi2S;
76+
77+
// histogram registry for bookkeeping
78+
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
79+
80+
static constexpr int nSpecies = 8;
81+
static constexpr int nParameters = 1;
82+
static const std::vector<std::string> particleNames;
83+
static const std::vector<int> particlePDGCodes;
84+
static const std::vector<std::string> parameterNames;
85+
static const int defaultParameters[nSpecies][nParameters];
86+
static constexpr std::string_view particleNamesConstExpr[] = {"EtaC1S", "JPsi", "ChiC0", "ChiC1",
87+
"hC", "ChiC2", "EtaC2S", "Psi2S"};
88+
89+
uint32_t enabledBits = 0;
90+
91+
Configurable<LabeledArray<int>> enableGeneratedInfo{"enableGeneratedInfo",
92+
{defaultParameters[0], nSpecies,
93+
nParameters, particleNames, parameterNames},
94+
"Fill generated particle histograms for each species. 0: no, 1: yes"};
95+
96+
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "p_{T} (GeV/c)"};
97+
ConfigurableAxis axisCentrality{"axisCentrality", {100, 0.0f, 100.0f}, "Centrality"};
98+
99+
ConfigurableAxis axisNVertices{"axisNVertices", {10, -0.5f, 9.5f}, "N(vertices)"};
100+
101+
std::vector<uint32_t> genEtaC1S;
102+
std::vector<uint32_t> genJPsi;
103+
std::vector<uint32_t> genChiC0;
104+
std::vector<uint32_t> genChiC1;
105+
std::vector<uint32_t> genHC;
106+
std::vector<uint32_t> genChiC2;
107+
std::vector<uint32_t> genEtaC2S;
108+
std::vector<uint32_t> genPsi2S;
109+
110+
// Preslice
111+
Preslice<aod::McParticles> mcParticlePerMcCollision = o2::aod::mcparticle::mcCollisionId;
112+
113+
void init(InitContext&)
114+
{
115+
// setup map for fast checking if enabled
116+
static_for<0, nSpecies - 1>([&](auto i) {
117+
constexpr int index = i.value;
118+
int f = enableGeneratedInfo->get(particleNames[index].c_str(), "Enable");
119+
if (f == 1) {
120+
bitset(enabledBits, index);
121+
}
122+
});
123+
124+
// Creation of histograms: MC generated
125+
for (Int_t i = 0; i < nSpecies; i++) {
126+
histos.add(Form("h2dGenerated%s", particleNames[i].data()), Form("h2dGenerated%s", particleNames[i].data()), kTH2D, {axisCentrality, axisPt});
127+
}
128+
129+
histos.add("h2dNVerticesVsCentrality", "h2dNVerticesVsCentrality", kTH2D, {axisCentrality, axisNVertices});
130+
131+
// reserve space for generated vectors if that process enabled
132+
auto hBinFinder = histos.get<TH2>(HIST("h2dGeneratedEtaC1S"));
133+
LOGF(info, "Binned generated processing enabled. Initialising with %i elements...", hBinFinder->GetNcells());
134+
genEtaC1S.resize(hBinFinder->GetNcells(), 0);
135+
genJPsi.resize(hBinFinder->GetNcells(), 0);
136+
genChiC0.resize(hBinFinder->GetNcells(), 0);
137+
genChiC1.resize(hBinFinder->GetNcells(), 0);
138+
genHC.resize(hBinFinder->GetNcells(), 0);
139+
genChiC2.resize(hBinFinder->GetNcells(), 0);
140+
genEtaC2S.resize(hBinFinder->GetNcells(), 0);
141+
genPsi2S.resize(hBinFinder->GetNcells(), 0);
142+
LOGF(info, "Binned generated processing: init done.");
143+
}
144+
145+
void processReconstructedSimulation(aod::McCollision const& /*mcCollision*/, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As, aod::FT0Mults>> const& collisions, aod::McParticles const& mcParticles)
146+
{
147+
// this process function also checks if a given collision was reconstructed and checks explicitly for splitting, etc
148+
// identify best-of collision
149+
int biggestNContribs = -1;
150+
float bestCentrality = 100.5;
151+
for (auto& collision : collisions) {
152+
if (biggestNContribs < collision.numContrib()) {
153+
biggestNContribs = collision.numContrib();
154+
bestCentrality = collision.centFT0M();
155+
}
156+
}
157+
histos.fill(HIST("h2dNVerticesVsCentrality"), bestCentrality, collisions.size());
158+
159+
for (auto& mcp : mcParticles) {
160+
if (TMath::Abs(mcp.y()) < 0.5 /* && mcp.isPhysicalPrimary()*/) {
161+
static_for<0, nSpecies - 1>([&](auto i) {
162+
constexpr int index = i.value;
163+
if (mcp.pdgCode() == particlePDGCodes[index] && bitcheck(enabledBits, index)) {
164+
histos.fill(HIST("h2dGenerated") + HIST(particleNamesConstExpr[index]), bestCentrality, mcp.pt());
165+
}
166+
});
167+
}
168+
}
169+
}
170+
171+
void processBinnedGenerated(soa::Join<aod::McCollisions, aod::McCollsExtra> const& mcCollisions, aod::McParticles const& mcParticlesEntireTable)
172+
{
173+
// set to zero
174+
std::fill(genEtaC1S.begin(), genEtaC1S.end(), 0);
175+
std::fill(genJPsi.begin(), genJPsi.end(), 0);
176+
std::fill(genChiC0.begin(), genChiC0.end(), 0);
177+
std::fill(genChiC1.begin(), genChiC1.end(), 0);
178+
std::fill(genHC.begin(), genHC.end(), 0);
179+
std::fill(genChiC2.begin(), genChiC2.end(), 0);
180+
std::fill(genEtaC2S.begin(), genEtaC2S.end(), 0);
181+
std::fill(genPsi2S.begin(), genPsi2S.end(), 0);
182+
183+
// this process function also checks if a given collision was reconstructed and checks explicitly for splitting, etc
184+
for (auto& mcCollision : mcCollisions) {
185+
const uint64_t mcCollIndex = mcCollision.globalIndex();
186+
187+
// use one of the generated histograms as the bin finder
188+
auto hBinFinder = histos.get<TH2>(HIST("h2dGeneratedEtaC1S"));
189+
190+
auto mcParticles = mcParticlesEntireTable.sliceBy(mcParticlePerMcCollision, mcCollIndex);
191+
for (auto& mcp : mcParticles) {
192+
if (TMath::Abs(mcp.y()) < 0.5 /* && mcp.isPhysicalPrimary()*/) {
193+
auto binNumber = hBinFinder->FindBin(mcCollision.bestCollisionCentFT0C(), mcp.pt()); // caution: pack
194+
if (mcp.pdgCode() == 441)
195+
genEtaC1S[binNumber]++;
196+
if (mcp.pdgCode() == 443)
197+
genJPsi[binNumber]++;
198+
if (mcp.pdgCode() == 10441)
199+
genChiC0[binNumber]++;
200+
if (mcp.pdgCode() == 20443)
201+
genChiC1[binNumber]++;
202+
if (mcp.pdgCode() == 10443)
203+
genHC[binNumber]++;
204+
if (mcp.pdgCode() == 445)
205+
genChiC2[binNumber]++;
206+
if (mcp.pdgCode() == 100441)
207+
genEtaC2S[binNumber]++;
208+
if (mcp.pdgCode() == 100443)
209+
genPsi2S[binNumber]++;
210+
}
211+
}
212+
}
213+
// at end of data frame
214+
// -> pack information from this DF into a generated histogram, once / DF
215+
geEtaC1S(genEtaC1S);
216+
geJPsi(genJPsi);
217+
geChiC0(genChiC0);
218+
geChiC1(genChiC1);
219+
geHC(genHC);
220+
geChiC2(genChiC2);
221+
geEtaC2S(genEtaC2S);
222+
gePsi2S(genPsi2S);
223+
}
224+
225+
PROCESS_SWITCH(generatedQuarkoniaMC, processReconstructedSimulation, "Produce reco-ed simulated information", true);
226+
PROCESS_SWITCH(generatedQuarkoniaMC, processBinnedGenerated, "Produce binned generated information", false);
227+
};
228+
229+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
230+
{
231+
return WorkflowSpec{
232+
adaptAnalysisTask<generatedQuarkoniaMC>(cfgc)};
233+
}
234+
235+
//__________________________________________________
236+
// do not over-populate general namespace, keep scope generatedQuarkoniaMC::
237+
const std::vector<std::string> generatedQuarkoniaMC::particleNames{"EtaC1S", "JPsi", "ChiC0", "ChiC1",
238+
"hC", "ChiC2", "EtaC2S", "Psi2S"};
239+
const std::vector<int> generatedQuarkoniaMC::particlePDGCodes{441, 443, 10441, 20443, 10443, 445, 100441, 100443};
240+
const std::vector<std::string> generatedQuarkoniaMC::parameterNames{"Enable"};
241+
242+
const int generatedQuarkoniaMC::defaultParameters[generatedQuarkoniaMC::nSpecies][generatedQuarkoniaMC::nParameters] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};

‎PWGDQ/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -107,4 +107,9 @@ o2physics_add_dpl_workflow(task-muon-mid-eff
107107
o2physics_add_dpl_workflow(task-fwd-track-pid
108108
SOURCES taskFwdTrackPid.cxx
109109
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore
110+
COMPONENT_NAME Analysis)
111+
112+
o2physics_add_dpl_workflow(quarkonia-to-hyperons
113+
SOURCES quarkoniaToHyperons.cxx
114+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore
110115
COMPONENT_NAME Analysis)

0 commit comments

Comments
 (0)