Skip to content

Commit ac265e6

Browse files
wang-yuanzhewang-yuanzhe
authored
[PWGLF] Add QA task for hyperhelium4sigma mc production (#11054)
1 parent 8c821c2 commit ac265e6

File tree

2 files changed

+343
-0
lines changed

2 files changed

+343
-0
lines changed

‎PWGLF/Tasks/Nuspex/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -149,6 +149,11 @@ o2physics_add_dpl_workflow(nuclei-from-hypertriton-map
149149
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
150150
COMPONENT_NAME Analysis)
151151

152+
o2physics_add_dpl_workflow(hyperhelium4sigma-analysis
153+
SOURCES hyperhelium4sigmaAnalysis.cxx
154+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
155+
COMPONENT_NAME Analysis)
156+
152157
if(FastJet_FOUND)
153158
o2physics_add_dpl_workflow(angular-correlations-in-jets
154159
SOURCES angularCorrelationsInJets.cxx

‎PWGLF/Tasks/Nuspex/hyperhelium4sigmaAnalysis.cxx‎

Lines changed: 338 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,338 @@
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 hyperhelium4sigmaAnalysis.cxx
13+
/// \brief Simple check for injected hyper-helium4sigma (H4S) in MC productions
14+
/// \author Yuanzhe Wang <yuanzhe.wang@cern.ch>
15+
16+
#include <vector>
17+
18+
#include "Framework/runDataProcessing.h"
19+
#include "Framework/AnalysisTask.h"
20+
#include "Framework/AnalysisDataModel.h"
21+
#include "Framework/ASoAHelpers.h"
22+
#include "ReconstructionDataFormats/Track.h"
23+
#include "Common/Core/RecoDecay.h"
24+
#include "Common/DataModel/EventSelection.h"
25+
#include "Common/DataModel/PIDResponse.h"
26+
#include "CommonConstants/PhysicsConstants.h"
27+
28+
using namespace o2;
29+
using namespace o2::framework;
30+
using namespace o2::framework::expressions;
31+
using std::array;
32+
using FullTracksExtIU = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::pidTPCFullPr, aod::pidTPCFullAl, aod::pidTPCFullTr, aod::pidTPCFullPi>;
33+
using MCLabeledTracksIU = soa::Join<FullTracksExtIU, aod::McTrackLabels>;
34+
35+
//-------------------------------Check the decay channel of H4S-------------------------------
36+
enum Channel {
37+
k2body = 0, // helium4, pion0
38+
k3body_p, // triton, proton, pion0
39+
k3body_n, // triton, neutron, pion+
40+
kNDecayChannel
41+
};
42+
43+
template <class TMCTrackTo, typename TMCParticle>
44+
Channel getDecayChannelH4S(TMCParticle const& particle)
45+
{
46+
if (std::abs(particle.pdgCode()) != o2::constants::physics::Pdg::kHyperHelium4Sigma) {
47+
return kNDecayChannel;
48+
}
49+
bool haveAlpha = false, haveTriton = false, haveProton = false, haveNeuteron = false;
50+
bool haveAntiAlpha = false, haveAntiTriton = false, haveAntiProton = false, haveAntiNeuteron = false;
51+
bool havePionPlus = false, havePionMinus = false, havePion0 = false;
52+
auto daughters = particle.template daughters_as<TMCTrackTo>();
53+
for (const auto& mcDaughter : particle.template daughters_as<TMCTrackTo>()) {
54+
if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kAlpha) {
55+
haveAlpha = true;
56+
}
57+
if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kAlpha) {
58+
haveAntiAlpha = true;
59+
}
60+
if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kTriton) {
61+
haveTriton = true;
62+
}
63+
if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kTriton) {
64+
haveAntiTriton = true;
65+
}
66+
if (mcDaughter.pdgCode() == PDG_t::kProton) {
67+
haveProton = true;
68+
}
69+
if (mcDaughter.pdgCode() == -PDG_t::kProton) {
70+
haveAntiProton = true;
71+
}
72+
if (mcDaughter.pdgCode() == PDG_t::kNeutron) {
73+
haveNeuteron = true;
74+
}
75+
if (mcDaughter.pdgCode() == -PDG_t::kNeutron) {
76+
haveAntiNeuteron = true;
77+
}
78+
if (mcDaughter.pdgCode() == PDG_t::kPiPlus) {
79+
havePionPlus = true;
80+
}
81+
if (mcDaughter.pdgCode() == -PDG_t::kPiPlus) {
82+
havePionMinus = true;
83+
}
84+
if (mcDaughter.pdgCode() == PDG_t::kPi0) {
85+
havePion0 = true;
86+
}
87+
}
88+
89+
if ((haveAlpha && havePion0) || (haveAntiAlpha && havePion0)) {
90+
return k2body;
91+
} else if ((haveTriton && haveProton && havePion0) || (haveAntiTriton && haveAntiProton && havePion0)) {
92+
return k3body_p;
93+
} else if ((haveTriton && haveNeuteron && havePionPlus) || (haveAntiTriton && haveAntiNeuteron && havePionMinus)) {
94+
return k3body_n;
95+
}
96+
97+
return kNDecayChannel;
98+
}
99+
//--------------------------------------------------------------
100+
101+
// check the performance of mcparticle
102+
struct Hyperhelium4sigmaAnalysis {
103+
// Basic checks
104+
HistogramRegistry registry{
105+
"registry",
106+
{
107+
{"hCollCounter", "hCollCounter", {HistType::kTH1F, {{2, 0.0f, 2.0f}}}},
108+
{"hMcCollCounter", "hMcCollCounter", {HistType::kTH1F, {{2, 0.0f, 2.0f}}}},
109+
},
110+
};
111+
112+
ConfigurableAxis ptBins{"ptBins", {200, 0.f, 10.f}, "Binning for #it{p}_{T} (GeV/#it{c})"};
113+
ConfigurableAxis ctBins{"ctBins", {100, 0.f, 25.f}, "Binning for c#it{t} (cm)"};
114+
ConfigurableAxis rigidityBins{"rigidityBins", {200, -10.f, 10.f}, "Binning for #it{p}_{T} (GeV/#it{c})"};
115+
ConfigurableAxis nsigmaBins{"nsigmaBins", {120, -6.f, 6.f}, "Binning for n sigma"};
116+
ConfigurableAxis invMassBins{"invMassBins", {100, 3.85, 4.15f}, "Binning for invariant mass (GeV/#it{c}^{2})"};
117+
118+
void init(InitContext&)
119+
{
120+
const AxisSpec ptAxis{ptBins, "#it{p}_{T} (GeV/#it{c})"};
121+
const AxisSpec ctAxis{ctBins, "c#it{t} (cm)"};
122+
const AxisSpec rigidityAxis{rigidityBins, "p/z (GeV/#it{c})"};
123+
const AxisSpec nsigmaAxis{nsigmaBins, "TPC n#sigma"};
124+
const AxisSpec invMassAxis{invMassBins, "Inv Mass (GeV/#it{c}^{2})"};
125+
126+
registry.get<TH1>(HIST("hCollCounter"))->GetXaxis()->SetBinLabel(1, "Reconstructed Collisions");
127+
registry.get<TH1>(HIST("hCollCounter"))->GetXaxis()->SetBinLabel(2, "Selected");
128+
registry.get<TH1>(HIST("hMcCollCounter"))->GetXaxis()->SetBinLabel(1, "MC Collisions");
129+
registry.get<TH1>(HIST("hMcCollCounter"))->GetXaxis()->SetBinLabel(2, "Reconstructed");
130+
131+
auto hGenHyperHelium4SigmaCounter = registry.add<TH1>("hGenHyperHelium4SigmaCounter", "", HistType::kTH1F, {{10, 0.f, 10.f}});
132+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(1, "H4S All");
133+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(2, "Matter");
134+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(3, "AntiMatter");
135+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(4, "#alpha + #pi^{0}");
136+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(5, "#bar{#alpha} + #pi^{0}");
137+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(6, "t + p + #pi^{0}");
138+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(7, "#bar{t} + #bar{p} + #pi^{0}");
139+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(8, "t + n + #pi^{+}");
140+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(9, "#bar{t} + #bar{n} + #pi^{+}");
141+
registry.get<TH1>(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(10, "Unexpected");
142+
143+
auto hEvtSelectedHyperHelium4SigmaCounter = registry.add<TH1>("hEvtSelectedHyperHelium4SigmaCounter", "", HistType::kTH1F, {{2, 0.f, 2.f}});
144+
registry.get<TH1>(HIST("hEvtSelectedHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(1, "Generated");
145+
registry.get<TH1>(HIST("hEvtSelectedHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(2, "Survived");
146+
147+
registry.add<TH1>("hGenHyperHelium4SigmaPt", "", HistType::kTH1F, {ptAxis});
148+
registry.add<TH1>("hGenHyperHelium4SigmaCt", "", HistType::kTH1F, {ctAxis});
149+
registry.add<TH1>("hMcRecoInvMass", "", HistType::kTH1F, {invMassAxis});
150+
151+
registry.add<TH2>("hDauHelium4TPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis});
152+
registry.add<TH2>("hDauTritonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis});
153+
registry.add<TH2>("hDauProtonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis});
154+
}
155+
156+
// Configurable<bool> eventSel8Cut{"eventSel8Cut", true, "flag to enable event sel8 selection"};
157+
Configurable<bool> mcEventCut{"mcEventCut", true, "flag to enable mc event selection: kIsTriggerTVX and kNoTimeFrameBorder"};
158+
Configurable<bool> eventPosZCut{"eventPosZCut", true, "flag to enable event posZ selection"};
159+
Configurable<float> maxPosZ{"maxPosZ", 10.f, "max pv posZ for event selection"};
160+
161+
Preslice<aod::McParticles> permcCollision = o2::aod::mcparticle::mcCollisionId;
162+
163+
std::vector<int64_t> mcPartIndices;
164+
template <typename TTrackTable>
165+
void setTrackIDForMC(aod::McParticles const& particlesMC, TTrackTable const& tracks)
166+
{
167+
mcPartIndices.clear();
168+
mcPartIndices.resize(particlesMC.size());
169+
std::fill(mcPartIndices.begin(), mcPartIndices.end(), -1);
170+
for (const auto& track : tracks) {
171+
if (track.has_mcParticle()) {
172+
auto mcparticle = track.template mcParticle_as<aod::McParticles>();
173+
if (mcPartIndices[mcparticle.globalIndex()] == -1) {
174+
mcPartIndices[mcparticle.globalIndex()] = track.globalIndex();
175+
} else {
176+
auto candTrack = tracks.rawIteratorAt(mcPartIndices[mcparticle.globalIndex()]);
177+
// Use the track which has innest information (also best quality?
178+
if (track.x() < candTrack.x()) {
179+
mcPartIndices[mcparticle.globalIndex()] = track.globalIndex();
180+
}
181+
}
182+
}
183+
}
184+
}
185+
186+
void process(aod::McCollisions const& mcCollisions, aod::McParticles const& particlesMC, o2::soa::Join<o2::aod::Collisions, o2::aod::McCollisionLabels, o2::aod::EvSels> const& collisions, MCLabeledTracksIU const& tracks)
187+
{
188+
setTrackIDForMC(particlesMC, tracks);
189+
std::vector<int64_t> selectedEvents(collisions.size());
190+
int nevts = 0;
191+
for (const auto& collision : collisions) {
192+
registry.fill(HIST("hCollCounter"), 0.5);
193+
if (mcEventCut && (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder))) {
194+
continue;
195+
}
196+
if (eventPosZCut && std::abs(collision.posZ()) > maxPosZ) { // 10cm
197+
continue;
198+
}
199+
registry.fill(HIST("hCollCounter"), 1.5);
200+
selectedEvents[nevts++] = collision.mcCollision_as<aod::McCollisions>().globalIndex();
201+
}
202+
selectedEvents.resize(nevts);
203+
204+
for (const auto& mcCollision : mcCollisions) {
205+
registry.fill(HIST("hMcCollCounter"), 0.5);
206+
const auto evtReconstructedAndSelected = std::find(selectedEvents.begin(), selectedEvents.end(), mcCollision.globalIndex()) != selectedEvents.end();
207+
if (evtReconstructedAndSelected) { // Check that the event is reconstructed and that the reconstructed events pass the selection
208+
registry.fill(HIST("hMcCollCounter"), 1.5);
209+
} else {
210+
// continue;
211+
}
212+
213+
const auto& dparticlesMC = particlesMC.sliceBy(permcCollision, mcCollision.globalIndex());
214+
215+
for (const auto& mcparticle : dparticlesMC) {
216+
217+
bool isMatter;
218+
if (mcparticle.pdgCode() == o2::constants::physics::Pdg::kHyperHelium4Sigma) {
219+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 1.5);
220+
isMatter = true;
221+
} else if (mcparticle.pdgCode() == -o2::constants::physics::Pdg::kHyperHelium4Sigma) {
222+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 2.5);
223+
isMatter = false;
224+
} else {
225+
continue;
226+
}
227+
228+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 0.5);
229+
registry.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 0.5);
230+
if (evtReconstructedAndSelected) {
231+
registry.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 1.5);
232+
}
233+
234+
double svPos[3] = {-999, -999, -999};
235+
double dauHelium4Mom[3] = {-999, -999, -999};
236+
double dauTritonMom[3] = {-999, -999, -999};
237+
double dauProtonMom[3] = {-999, -999, -999};
238+
double dauNeuteronMom[3] = {-999, -999, -999};
239+
double dauChargedPionMom[3] = {-999, -999, -999};
240+
double dauPion0Mom[3] = {-999, -999, -999};
241+
auto dChannel = getDecayChannelH4S<aod::McParticles>(mcparticle);
242+
if (dChannel == kNDecayChannel) {
243+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 9.5);
244+
continue;
245+
}
246+
for (const auto& mcparticleDaughter : mcparticle.daughters_as<aod::McParticles>()) {
247+
if (std::abs(mcparticleDaughter.pdgCode()) == o2::constants::physics::Pdg::kAlpha) {
248+
dauHelium4Mom[0] = mcparticleDaughter.px();
249+
dauHelium4Mom[1] = mcparticleDaughter.py();
250+
dauHelium4Mom[2] = mcparticleDaughter.pz();
251+
252+
// get SV position for 2body decay
253+
svPos[0] = mcparticleDaughter.vx();
254+
svPos[1] = mcparticleDaughter.vy();
255+
svPos[2] = mcparticleDaughter.vz();
256+
257+
if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) {
258+
auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]);
259+
registry.fill(HIST("hDauHelium4TPCNSigma"), track.p() * track.sign(), track.tpcNSigmaAl());
260+
}
261+
} else if (std::abs(mcparticleDaughter.pdgCode()) == o2::constants::physics::Pdg::kTriton) {
262+
dauTritonMom[0] = mcparticleDaughter.px();
263+
dauTritonMom[1] = mcparticleDaughter.py();
264+
dauTritonMom[2] = mcparticleDaughter.pz();
265+
266+
// get SV position for 3body decay
267+
svPos[0] = mcparticleDaughter.vx();
268+
svPos[1] = mcparticleDaughter.vy();
269+
svPos[2] = mcparticleDaughter.vz();
270+
271+
if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) {
272+
auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]);
273+
registry.fill(HIST("hDauTritonTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaTr());
274+
}
275+
} else if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kProton) {
276+
dauProtonMom[0] = mcparticleDaughter.px();
277+
dauProtonMom[1] = mcparticleDaughter.py();
278+
dauProtonMom[2] = mcparticleDaughter.pz();
279+
280+
if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) {
281+
auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]);
282+
registry.fill(HIST("hDauProtonTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaPr());
283+
}
284+
} else if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kNeutron) {
285+
dauNeuteronMom[0] = mcparticleDaughter.px();
286+
dauNeuteronMom[1] = mcparticleDaughter.py();
287+
dauNeuteronMom[2] = mcparticleDaughter.pz();
288+
} else if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kPiPlus) {
289+
dauChargedPionMom[0] = mcparticleDaughter.px();
290+
dauChargedPionMom[1] = mcparticleDaughter.py();
291+
dauChargedPionMom[2] = mcparticleDaughter.pz();
292+
} else if (mcparticleDaughter.pdgCode() == PDG_t::kPi0) {
293+
dauPion0Mom[0] = mcparticleDaughter.px();
294+
dauPion0Mom[1] = mcparticleDaughter.py();
295+
dauPion0Mom[2] = mcparticleDaughter.pz();
296+
}
297+
}
298+
299+
registry.fill(HIST("hGenHyperHelium4SigmaPt"), mcparticle.pt());
300+
double ct = RecoDecay::sqrtSumOfSquares(svPos[0] - mcparticle.vx(), svPos[1] - mcparticle.vy(), svPos[2] - mcparticle.vz()) * o2::constants::physics::MassHyperHelium4Sigma / mcparticle.p();
301+
registry.fill(HIST("hGenHyperHelium4SigmaCt"), ct);
302+
303+
if (dChannel == k2body) {
304+
if (isMatter) {
305+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 3.5);
306+
} else {
307+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 4.5);
308+
}
309+
double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauHelium4Mom[0], dauHelium4Mom[1], dauHelium4Mom[2]}, std::array{dauPion0Mom[0], dauPion0Mom[1], dauPion0Mom[2]}}, std::array{o2::constants::physics::MassAlpha, o2::constants::physics::MassPi0});
310+
registry.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass);
311+
} else if (dChannel == k3body_p) {
312+
if (isMatter) {
313+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 5.5);
314+
} else {
315+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 6.5);
316+
}
317+
double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauTritonMom[0], dauTritonMom[1], dauTritonMom[2]}, std::array{dauProtonMom[0], dauProtonMom[1], dauProtonMom[2]}, std::array{dauPion0Mom[0], dauPion0Mom[1], dauPion0Mom[2]}}, std::array{o2::constants::physics::MassTriton, o2::constants::physics::MassProton, o2::constants::physics::MassPi0});
318+
registry.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass);
319+
} else if (dChannel == k3body_n) {
320+
if (isMatter) {
321+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 7.5);
322+
} else {
323+
registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 8.5);
324+
}
325+
double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauTritonMom[0], dauTritonMom[1], dauTritonMom[2]}, std::array{dauNeuteronMom[0], dauNeuteronMom[1], dauNeuteronMom[2]}, std::array{dauChargedPionMom[0], dauChargedPionMom[1], dauChargedPionMom[2]}}, std::array{o2::constants::physics::MassTriton, o2::constants::physics::MassNeutron, o2::constants::physics::MassPionCharged});
326+
registry.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass);
327+
}
328+
}
329+
}
330+
}
331+
};
332+
333+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
334+
{
335+
return WorkflowSpec{
336+
adaptAnalysisTask<Hyperhelium4sigmaAnalysis>(cfgc),
337+
};
338+
}

0 commit comments

Comments
 (0)