Skip to content

Commit fc7b1e9

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

File tree

2 files changed

+201
-1
lines changed

2 files changed

+201
-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: 195 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,195 @@
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+
// using BCs = soa::Join<aod::BCsWithTimestamps, aod::BcSels, aod::Run3MatchedToBCSparse>;
53+
54+
Configurable<float> minB{"minB", 0.0f, "min impact parameter"};
55+
Configurable<float> maxB{"maxB", 20.0f, "max impact parameter"};
56+
O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
57+
O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
58+
O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.1f, "Minimal pT for tracks")
59+
O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 1000.0f, "Maximal pT for tracks")
60+
O2_DEFINE_CONFIGURABLE(cfgCutDCAxy, float, 0.2f, "DCAxy cut for tracks")
61+
O2_DEFINE_CONFIGURABLE(cfgDcaxy, bool, true, "choose dcaxy")
62+
63+
ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""};
64+
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"};
65+
// Connect to ccdb
66+
Service<ccdb::BasicCCDBManager> ccdb;
67+
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
68+
69+
double epsilon = 1e-6;
70+
71+
using McParticles = soa::Join<aod::UDMcParticles, aod::UDMcTrackLabels>;
72+
73+
void init(InitContext&)
74+
{
75+
ccdb->setURL(ccdbUrl.value);
76+
ccdb->setCaching(true);
77+
auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
78+
ccdb->setCreatedNotAfter(now);
79+
80+
const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
81+
const AxisSpec axisEta{20, -1., 1., "#eta"};
82+
const AxisSpec axisCounter{1, 0, +1, ""};
83+
// QA histograms
84+
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {{5, 0, 5}});
85+
histos.add<TH1>("RecoProcessEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}});
86+
histos.add<TH1>("hImpactParameter", "hImpactParameter", HistType::kTH1D, {axisB});
87+
88+
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
89+
histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
90+
histos.add<TH1>("hPtReco", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
91+
histos.add<TH3>("hEtaPtVtxzMCReco", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
92+
}
93+
94+
template <typename TCollision>
95+
bool eventSelected(TCollision collision)
96+
{
97+
return true;
98+
}
99+
100+
template <typename TTrack>
101+
bool trackSelected(TTrack track)
102+
{
103+
auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
104+
double pt = RecoDecay::pt(momentum);
105+
if (pt < cfgPtCutMin || pt > cfgPtCutMax) {
106+
return false;
107+
}
108+
double dcaLimit = 0.0105 + 0.035 / std::pow(pt, 1.1);
109+
if (cfgDcaxy && !(std::fabs(track.dcaXY()) < dcaLimit)) {
110+
return false;
111+
}
112+
return true;
113+
}
114+
115+
void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, McParticles const& mcParticles)
116+
{
117+
// if (bcs.size() == 0) {
118+
// return;
119+
// }
120+
histos.fill(HIST("mcEventCounter"), 0.5);
121+
float imp = mcCollision.impactParameter();
122+
float vtxz = mcCollision.posZ();
123+
124+
if (imp >= minB && imp <= maxB) {
125+
// event within range
126+
histos.fill(HIST("hImpactParameter"), imp);
127+
128+
for (auto const& mcParticle : mcParticles) {
129+
auto momentum = std::array<double, 3>{mcParticle.px(), mcParticle.py(), mcParticle.pz()};
130+
double pt = RecoDecay::pt(momentum);
131+
double phi = RecoDecay::phi(momentum);
132+
double eta = RecoDecay::eta(momentum);
133+
// focus on bulk: e, mu, pi, k, p
134+
int pdgCode = std::abs(mcParticle.pdgCode());
135+
if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
136+
continue;
137+
138+
if (!mcParticle.isPhysicalPrimary())
139+
continue;
140+
// if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
141+
// continue;
142+
143+
histos.fill(HIST("hPtMCGen"), pt);
144+
histos.fill(HIST("hEtaPtVtxzMCGen"), eta, pt, vtxz);
145+
}
146+
}
147+
}
148+
PROCESS_SWITCH(FlowMcUpc, processMCTrue, "process pure simulation information", true);
149+
150+
using MCRecoTracks = soa::Join<aod::UDTracks, aod::UDTracksPID, aod::UDTracksExtra, aod::UDTracksFlags, aod::UDTracksDCA, aod::UDMcTrackLabels>;
151+
using MCRecoCollisions = soa::Join<aod::UDCollisions, aod::UDCollisionsSels, aod::UDMcCollsLabels>;
152+
153+
void processReco(MCRecoCollisions::iterator const& collision, MCRecoTracks const& tracks)
154+
{
155+
histos.fill(HIST("RecoProcessEventCounter"), 0.5);
156+
if (!eventSelected(collision))
157+
return;
158+
histos.fill(HIST("RecoProcessEventCounter"), 1.5);
159+
if (!collision.has_udMcCollision())
160+
return;
161+
histos.fill(HIST("RecoProcessEventCounter"), 2.5);
162+
if (tracks.size() < 1)
163+
return;
164+
histos.fill(HIST("RecoProcessEventCounter"), 3.5);
165+
166+
float vtxz = collision.posZ();
167+
168+
for (const auto& track : tracks) {
169+
auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
170+
double pt = RecoDecay::pt(momentum);
171+
double phi = RecoDecay::phi(momentum);
172+
double eta = RecoDecay::eta(momentum);
173+
if (!trackSelected(track) || (!track.has_udMcParticle()))
174+
continue;
175+
auto mcParticle = track.udMcParticle();
176+
int pdgCode = std::abs(mcParticle.pdgCode());
177+
if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
178+
continue;
179+
// if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
180+
// continue;
181+
if (!mcParticle.isPhysicalPrimary())
182+
continue;
183+
184+
histos.fill(HIST("hPtReco"), pt);
185+
histos.fill(HIST("hEtaPtVtxzMCReco"), eta, pt, vtxz);
186+
}
187+
}
188+
PROCESS_SWITCH(FlowMcUpc, processReco, "process reconstructed information", true);
189+
};
190+
191+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
192+
{
193+
return WorkflowSpec{
194+
adaptAnalysisTask<FlowMcUpc>(cfgc)};
195+
}

0 commit comments

Comments
 (0)