Skip to content

Commit 63792d9

Browse files
fchinufchinu
committed
Add checkMcPvContr task
1 parent b179555 commit 63792d9

2 files changed

Lines changed: 177 additions & 0 deletions

File tree

‎PWGHF/Tasks/CMakeLists.txt‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,16 @@ o2physics_add_dpl_workflow(task-pid-studies
4949
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils O2Physics::AnalysisCCDB
5050
COMPONENT_NAME Analysis)
5151

52+
o2physics_add_dpl_workflow(check-mc-pv-contr
53+
SOURCES checkMcPvContr.cxx
54+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils
55+
COMPONENT_NAME Analysis)
56+
57+
o2physics_add_dpl_workflow(check-aod
58+
SOURCES checkAod.cxx
59+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils O2Physics::AnalysisCCDB
60+
COMPONENT_NAME Analysis)
61+
5262
# o2physics_add_dpl_workflow(task-sel-optimisation
5363
# SOURCES taskSelOptimisation.cxx
5464
# PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore

‎PWGHF/Tasks/checkMcPvContr.cxx‎

Lines changed: 167 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,167 @@
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 taskMultiplicityEstimatorCorrelation.cxx
13+
/// \brief Task for correlating the multiplicity estimator with generated dN/deta
14+
///
15+
/// \author Fabrizio Chinu <fabrizio.chinu@cern.ch>, Università and INFN Torino
16+
17+
#include "Common/DataModel/Centrality.h"
18+
#include "Common/DataModel/EventSelection.h"
19+
#include "Common/DataModel/Multiplicity.h"
20+
21+
#include <Framework/ASoA.h>
22+
#include <Framework/AnalysisDataModel.h>
23+
#include <Framework/AnalysisTask.h>
24+
#include <Framework/Configurable.h>
25+
#include <Framework/HistogramRegistry.h>
26+
#include <Framework/HistogramSpec.h>
27+
#include <Framework/InitContext.h>
28+
#include <Framework/StaticFor.h>
29+
#include <Framework/runDataProcessing.h>
30+
31+
#include <TH2.h>
32+
#include <TPDGCode.h>
33+
34+
#include <array>
35+
#include <cstdint>
36+
#include <cstdlib>
37+
#include <memory>
38+
#include <string>
39+
#include <string_view>
40+
#include <vector>
41+
42+
using namespace o2;
43+
using namespace o2::framework;
44+
using namespace o2::framework::expressions;
45+
46+
namespace o2::aod
47+
{
48+
namespace check_mc_pv_contr
49+
{
50+
// V0s
51+
DECLARE_SOA_COLUMN(IdCollGen, idCollGen, int); //! Generated collision index
52+
DECLARE_SOA_COLUMN(ImpParGen, impParGen, float); //! Generated impact parameter
53+
DECLARE_SOA_COLUMN(TimeGen, timeGen, float); //! Generated collision time
54+
DECLARE_SOA_COLUMN(TimeRec, timeRec, float); //! Reconstructed collision time
55+
DECLARE_SOA_COLUMN(NCharm, nCharm, int); //! Number of charm quarks in the collision
56+
DECLARE_SOA_COLUMN(NCharmFromInj, nCharmFromInj, int); //! Number of charm quarks from injected events
57+
DECLARE_SOA_COLUMN(NPVContributors, nPVContributors, int); //! Number of contributors to the PV
58+
DECLARE_SOA_COLUMN(Centrality, centrality, int); //! Centrality FT0C
59+
DECLARE_SOA_COLUMN(BC, Bc, int); //! Bunch crossing
60+
} // namespace check_mc_pv_contr
61+
62+
DECLARE_SOA_TABLE(CheckInj, "AOD", "CHECKINJ", //! Table with PID information
63+
check_mc_pv_contr::IdCollGen,
64+
check_mc_pv_contr::ImpParGen,
65+
check_mc_pv_contr::TimeGen,
66+
check_mc_pv_contr::TimeRec,
67+
check_mc_pv_contr::NCharm,
68+
check_mc_pv_contr::NCharmFromInj,
69+
check_mc_pv_contr::NPVContributors,
70+
check_mc_pv_contr::Centrality,
71+
check_mc_pv_contr::BC);
72+
} // namespace o2::aod
73+
74+
struct checkMcPvContr {
75+
76+
Produces<o2::aod::CheckInj> checkInj;
77+
78+
using TrackWLabels = soa::Join<aod::Tracks, aod::McTrackLabels>;
79+
using CollisionWLabels = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::CentFT0Cs, aod::CentNTPVs>;
80+
81+
PresliceUnsorted<CollisionWLabels> colPerMcCollision = aod::mccollisionlabel::mcCollisionId;
82+
Preslice<TrackWLabels> tracksPerCollision = aod::track::collisionId;
83+
84+
HistogramRegistry registry{"registry", {}};
85+
std::shared_ptr<TH1> hCharmPerCollImpPar;
86+
87+
void init(InitContext&)
88+
{
89+
registry.add("hCharmImpPar", ";Impact parameter (fm);Charm counts", {HistType::kTH1F, {{200, 0, 20}}});
90+
registry.add("hCollImpPar", ";Impact parameter (fm);Counts", {HistType::kTH1F, {{200, 0, 20}}});
91+
hCharmPerCollImpPar = registry.add<TH1>("hCharmPerCollImpPar", ";Impact parameter (fm);Charm counts per collision", {HistType::kTH1F, {{200, 0, 20}}});
92+
93+
registry.add("hDeltaX", ";#DeltaX (cm);Counts", {HistType::kTH1F, {{200, 0.01, 0.01}}});
94+
registry.add("hDeltaY", ";#DeltaY (cm);Counts", {HistType::kTH1F, {{200, 0.01, 0.01}}});
95+
registry.add("hDeltaZ", ";#DeltaZ (cm);Counts", {HistType::kTH1F, {{200, -0.01, 0.01}}});
96+
}
97+
98+
void process(CollisionWLabels const& collisions,
99+
TrackWLabels const& tracks,
100+
aod::McParticles const& mcParticles,
101+
aod::McCollisions const& mcCollisions)
102+
{
103+
int splitColls{0};
104+
for (const auto& mcColl : mcCollisions) {
105+
const auto collSlice = collisions.sliceBy(colPerMcCollision, mcColl.globalIndex());
106+
int64_t idxCollMostPVContrib{-1};
107+
int nPVContribMost{0};
108+
// First we find the collision with most PV contributors
109+
for (const auto& collision : collSlice) {
110+
if (collision.centFT0C() < 20.f) {
111+
if (collision.numContrib() > nPVContribMost) {
112+
nPVContribMost = collision.numContrib();
113+
idxCollMostPVContrib = collision.globalIndex();
114+
}
115+
}
116+
}
117+
118+
// Then we fill the histogram with the distances of the collisions
119+
for (const auto& collision : collSlice) {
120+
const auto collTracks = tracks.sliceBy(tracksPerCollision, collision.globalIndex());
121+
std::vector<int> charmIds{};
122+
int fromSignalEv{0};
123+
for (const auto& track : collTracks) {
124+
if (track.has_mcParticle()) {
125+
auto mcPart = track.mcParticle_as<aod::McParticles>();
126+
for (const auto& mother : mcPart.mothers_as<aod::McParticles>()) {
127+
if (std::abs(mother.pdgCode()) == 4) { // (anti)charm quark
128+
registry.fill(HIST("hDeltaX"), collision.posX() - collisions.rawIteratorAt(idxCollMostPVContrib).posX());
129+
registry.fill(HIST("hDeltaY"), collision.posY() - collisions.rawIteratorAt(idxCollMostPVContrib).posY());
130+
registry.fill(HIST("hDeltaZ"), collision.posZ() - collisions.rawIteratorAt(idxCollMostPVContrib).posZ());
131+
if (std::find(charmIds.begin(), charmIds.end(), mother.globalIndex()) == charmIds.end()) {
132+
charmIds.push_back(mother.globalIndex());
133+
fromSignalEv += static_cast<int>(!mother.fromBackgroundEvent());
134+
}
135+
break;
136+
}
137+
}
138+
}
139+
}
140+
checkInj(mcColl.globalIndex(), mcColl.impactParameter(), mcColl.t(), collision.collisionTime(), charmIds.size(), fromSignalEv, collision.numContrib(), collision.centFT0C(), collision.bcId());
141+
}
142+
}
143+
for (const auto& mcColl : mcCollisions) {
144+
registry.fill(HIST("hCollImpPar"), mcColl.impactParameter());
145+
}
146+
int count_bkg{0}, count_sgn{0};
147+
for (const auto& mcPart : mcParticles) {
148+
if (std::abs(mcPart.pdgCode()) == 4) { // (anti)charm quark
149+
if (!mcPart.fromBackgroundEvent()) {
150+
auto mcCollision = mcPart.mcCollision_as<aod::McCollisions>();
151+
registry.fill(HIST("hCharmImpPar"), mcCollision.impactParameter());
152+
count_sgn++;
153+
} else {
154+
count_bkg++;
155+
}
156+
}
157+
}
158+
hCharmPerCollImpPar->Divide(registry.get<TH1>(HIST("hCharmImpPar")).get(), registry.get<TH1>(HIST("hCollImpPar")).get(), 1, 1);
159+
LOG(info) << "Number of bkgev particles: " << count_bkg;
160+
LOG(info) << "Number of sngev particles: " << count_sgn;
161+
}
162+
};
163+
164+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
165+
{
166+
return WorkflowSpec{adaptAnalysisTask<checkMcPvContr>(cfgc)};
167+
}

0 commit comments

Comments
 (0)