Skip to content

Commit 7290de7

Browse files
miedema-11miedema-11
committed
MC analysis for flow
1 parent 913605e commit 7290de7

File tree

2 files changed

+199
-1
lines changed

2 files changed

+199
-1
lines changed

‎PWGUD/Tasks/CMakeLists.txt‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -254,6 +254,11 @@ o2physics_add_dpl_workflow(flow-correlations-upc
254254
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGCFCore
255255
COMPONENT_NAME Analysis)
256256

257+
o2physics_add_dpl_workflow(flow-mc-upc
258+
SOURCES flowMcUpc.cxx
259+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore
260+
COMPONENT_NAME Analysis)
261+
257262
o2physics_add_dpl_workflow(analysis-mc-dpm-jet-sg-v3
258263
SOURCES analysisMCDPMJetSGv3.cxx
259264
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
@@ -272,4 +277,4 @@ o2physics_add_dpl_workflow(sg-exclusive-jpsi-midrapidity
272277
o2physics_add_dpl_workflow(fitbit-mapping
273278
SOURCES upcTestFITBitMapping.cxx
274279
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
275-
COMPONENT_NAME Analysis)
280+
COMPONENT_NAME Analysis)

‎PWGUD/Tasks/flowMcUpc.cxx‎

Lines changed: 193 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,193 @@
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+
/// \file flowMcUpc.cxx
13+
/// \author Zhiyong Lu (zhiyong.lu@cern.ch), Yongxi Du (yongxi.du@cern.ch)
14+
/// \since Apr/2/2026
15+
/// \brief flow efficiency analysis on UPC MC
16+
17+
#include "PWGUD/Core/SGSelector.h"
18+
#include "PWGUD/DataModel/UDTables.h"
19+
20+
#include "Common/Core/RecoDecay.h"
21+
#include "Common/Core/TrackSelection.h"
22+
#include "Common/Core/TrackSelectionDefaults.h"
23+
#include "Common/Core/trackUtilities.h"
24+
#include "Common/DataModel/TrackSelectionTables.h"
25+
26+
#include "Framework/ASoAHelpers.h"
27+
#include "Framework/AnalysisDataModel.h"
28+
#include "Framework/AnalysisTask.h"
29+
#include "Framework/HistogramRegistry.h"
30+
#include "Framework/RunningWorkflowInfo.h"
31+
#include "Framework/runDataProcessing.h"
32+
#include "ReconstructionDataFormats/Track.h"
33+
#include <CCDB/BasicCCDBManager.h>
34+
35+
#include <TF1.h>
36+
#include <TPDGCode.h>
37+
#include <TProfile.h>
38+
#include <TRandom3.h>
39+
40+
#include <string>
41+
#include <vector>
42+
43+
using namespace o2;
44+
using namespace o2::framework;
45+
using namespace o2::framework::expressions;
46+
47+
#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
48+
49+
struct FlowMcUpc {
50+
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
51+
52+
Configurable<float> minB{"minB", 0.0f, "min impact parameter"};
53+
Configurable<float> maxB{"maxB", 20.0f, "max impact parameter"};
54+
O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
55+
O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
56+
O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.1f, "Minimal pT for tracks")
57+
O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 1000.0f, "Maximal pT for tracks")
58+
O2_DEFINE_CONFIGURABLE(cfgCutDCAxy, float, 0.2f, "DCAxy cut for tracks")
59+
O2_DEFINE_CONFIGURABLE(cfgDcaxy, bool, true, "choose dcaxy")
60+
61+
ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""};
62+
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}, "pt axis"};
63+
// Connect to ccdb
64+
Service<ccdb::BasicCCDBManager> ccdb;
65+
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
66+
67+
double epsilon = 1e-6;
68+
69+
using McParticles = soa::Join<aod::UDMcParticles, aod::UDMcTrackLabels>;
70+
71+
void init(InitContext&)
72+
{
73+
ccdb->setURL(ccdbUrl.value);
74+
ccdb->setCaching(true);
75+
auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
76+
ccdb->setCreatedNotAfter(now);
77+
78+
const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
79+
const AxisSpec axisEta{20, -1., 1., "#eta"};
80+
const AxisSpec axisCounter{1, 0, +1, ""};
81+
// QA histograms
82+
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {{5, 0, 5}});
83+
histos.add<TH1>("RecoProcessEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}});
84+
histos.add<TH1>("hImpactParameter", "hImpactParameter", HistType::kTH1D, {axisB});
85+
86+
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
87+
histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
88+
histos.add<TH1>("hPtReco", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
89+
histos.add<TH3>("hEtaPtVtxzMCReco", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
90+
}
91+
92+
// template <typename TCollision>
93+
// bool eventSelected(TCollision collision)
94+
// {
95+
// return true;
96+
// }
97+
98+
template <typename TTrack>
99+
bool trackSelected(TTrack track)
100+
{
101+
auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
102+
double pt = RecoDecay::pt(momentum);
103+
if (pt < cfgPtCutMin || pt > cfgPtCutMax) {
104+
return false;
105+
}
106+
double dcaLimit = 0.0105 + 0.035 / std::pow(pt, 1.1);
107+
if (cfgDcaxy && !(std::fabs(track.dcaXY()) < dcaLimit)) {
108+
return false;
109+
}
110+
return true;
111+
}
112+
113+
void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, McParticles const& mcParticles, aod::BCs const& bcs)
114+
{
115+
if (bcs.size() == 0) {
116+
return;
117+
}
118+
histos.fill(HIST("mcEventCounter"), 0.5);
119+
float imp = mcCollision.impactParameter();
120+
float vtxz = mcCollision.posZ();
121+
122+
if (imp >= minB && imp <= maxB) {
123+
// event within range
124+
histos.fill(HIST("hImpactParameter"), imp);
125+
126+
for (auto const& mcParticle : mcParticles) {
127+
auto momentum = std::array<double, 3>{mcParticle.px(), mcParticle.py(), mcParticle.pz()};
128+
double pt = RecoDecay::pt(momentum);
129+
double phi = RecoDecay::phi(momentum);
130+
double eta = RecoDecay::eta(momentum);
131+
// focus on bulk: e, mu, pi, k, p
132+
int pdgCode = std::abs(mcParticle.pdgCode());
133+
if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
134+
continue;
135+
136+
if (!mcParticle.isPhysicalPrimary())
137+
continue;
138+
// if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
139+
// continue;
140+
141+
histos.fill(HIST("hPtMCGen"), pt);
142+
histos.fill(HIST("hEtaPtVtxzMCGen"), eta, pt, vtxz);
143+
}
144+
}
145+
}
146+
PROCESS_SWITCH(FlowMcUpc, processMCTrue, "process pure simulation information", true);
147+
148+
using MCRecoTracks = soa::Join<aod::UDTracks, aod::UDTracksPID, aod::UDTracksExtra, aod::UDTracksFlags, aod::UDTracksDCA, aod::UDMcTrackLabels>;
149+
using MCRecoCollisions = soa::Join<aod::UDCollisions, aod::SGCollisions, aod::UDCollisionSelExtras, aod::UDCollisionsSels, aod::UDZdcsReduced, aod::UDMcCollsLabels>;
150+
151+
void processReco(MCRecoCollisions::iterator const& collision, MCRecoTracks const& tracks)
152+
{
153+
histos.fill(HIST("RecoProcessEventCounter"), 0.5);
154+
// if (!eventSelected(collision))
155+
// return;
156+
histos.fill(HIST("RecoProcessEventCounter"), 1.5);
157+
if (!collision.has_udMcCollision())
158+
return;
159+
histos.fill(HIST("RecoProcessEventCounter"), 2.5);
160+
if (tracks.size() < 1)
161+
return;
162+
histos.fill(HIST("RecoProcessEventCounter"), 3.5);
163+
164+
float vtxz = collision.posZ();
165+
166+
for (const auto& track : tracks) {
167+
auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
168+
double pt = RecoDecay::pt(momentum);
169+
double phi = RecoDecay::phi(momentum);
170+
double eta = RecoDecay::eta(momentum);
171+
if (!trackSelected(track) || (!track.has_udMcParticle()))
172+
continue;
173+
auto mcParticle = track.udMcParticle();
174+
int pdgCode = std::abs(mcParticle.pdgCode());
175+
if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
176+
continue;
177+
// if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
178+
// continue;
179+
if (!mcParticle.isPhysicalPrimary())
180+
continue;
181+
182+
histos.fill(HIST("hPtReco"), pt);
183+
histos.fill(HIST("hEtaPtVtxzMCReco"), eta, pt, vtxz);
184+
}
185+
}
186+
PROCESS_SWITCH(FlowMcUpc, processReco, "process reconstructed information", true);
187+
};
188+
189+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
190+
{
191+
return WorkflowSpec{
192+
adaptAnalysisTask<FlowMcUpc>(cfgc)};
193+
}

0 commit comments

Comments
 (0)