Skip to content

Commit 8924407

Browse files
fchinufchinu
committed
Add more output information
1 parent 96f9ae4 commit 8924407

2 files changed

Lines changed: 100 additions & 57 deletions

File tree

‎PWGHF/Tasks/CMakeLists.txt‎

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -49,16 +49,11 @@ 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
52+
o2physics_add_dpl_workflow(task-mc-injection
53+
SOURCES taskMcInjection.cxx
5454
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils
5555
COMPONENT_NAME Analysis)
5656

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-
6257
# o2physics_add_dpl_workflow(task-sel-optimisation
6358
# SOURCES taskSelOptimisation.cxx
6459
# PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore

PWGHF/Tasks/checkMcPvContr.cxx renamed to PWGHF/Tasks/taskMcInjection.cxx

Lines changed: 98 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -9,8 +9,8 @@
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
1111

12-
/// \file taskMultiplicityEstimatorCorrelation.cxx
13-
/// \brief Task for correlating the multiplicity estimator with generated dN/deta
12+
/// \file taskMcInjection.cxx
13+
/// \brief Task for checking injected events in Pb-Pb MC productions
1414
///
1515
/// \author Fabrizio Chinu <fabrizio.chinu@cern.ch>, Università and INFN Torino
1616

@@ -47,7 +47,7 @@ namespace o2::aod
4747
{
4848
namespace check_mc_pv_contr
4949
{
50-
// V0s
50+
// Collisions
5151
DECLARE_SOA_COLUMN(IdCollGen, idCollGen, int); //! Generated collision index
5252
DECLARE_SOA_COLUMN(ImpParGen, impParGen, float); //! Generated impact parameter
5353
DECLARE_SOA_COLUMN(XCollGen, xCollGen, float); //! Generated x coordinate of the collision
@@ -63,6 +63,11 @@ DECLARE_SOA_COLUMN(XCollRec, xCollRec, float); //! Reconstructed x c
6363
DECLARE_SOA_COLUMN(YCollRec, yCollRec, float); //! Reconstructed y coordinate of the collision
6464
DECLARE_SOA_COLUMN(ZCollRec, zCollRec, float); //! Reconstructed z coordinate of the collision
6565
DECLARE_SOA_COLUMN(BC, Bc, int); //! Bunch crossing
66+
// Tracks
67+
DECLARE_SOA_COLUMN(VX, vx, float); // x coordinate of the track production vertex
68+
DECLARE_SOA_COLUMN(VY, vy, float); // y coordinate of the track production vertex
69+
DECLARE_SOA_COLUMN(VZ, vz, float); // z coordinate of the track production vertex
70+
DECLARE_SOA_COLUMN(IsFromSignal, isFromSignal, bool); // Whether the track is from the signal event
6671
} // namespace check_mc_pv_contr
6772

6873
DECLARE_SOA_TABLE(CheckInj, "AOD", "CHECKINJ", //! Table with PID information
@@ -72,29 +77,42 @@ DECLARE_SOA_TABLE(CheckInj, "AOD", "CHECKINJ", //! Table with PID information
7277
check_mc_pv_contr::YCollGen,
7378
check_mc_pv_contr::ZCollGen,
7479
check_mc_pv_contr::TimeGen,
80+
check_mc_pv_contr::XCollRec,
81+
check_mc_pv_contr::YCollRec,
82+
check_mc_pv_contr::ZCollRec,
7583
check_mc_pv_contr::TimeRec,
7684
check_mc_pv_contr::NCharm,
7785
check_mc_pv_contr::NCharmFromInj,
7886
check_mc_pv_contr::NPVContributors,
7987
check_mc_pv_contr::Centrality,
80-
check_mc_pv_contr::XCollRec,
81-
check_mc_pv_contr::YCollRec,
82-
check_mc_pv_contr::ZCollRec,
8388
check_mc_pv_contr::BC);
89+
90+
DECLARE_SOA_TABLE(TracksInjection, "AOD", "TRKINJ", //! Table with MC labels for tracks
91+
check_mc_pv_contr::IdCollGen,
92+
check_mc_pv_contr::VX,
93+
check_mc_pv_contr::VY,
94+
check_mc_pv_contr::VZ,
95+
check_mc_pv_contr::IsFromSignal);
8496
} // namespace o2::aod
8597

86-
struct checkMcPvContr {
98+
struct taskMcInjection {
8799

88100
Produces<o2::aod::CheckInj> checkInj;
101+
Produces<o2::aod::TracksInjection> tracksInj;
89102

90103
using TrackWLabels = soa::Join<aod::Tracks, aod::McTrackLabels>;
91104
using CollisionWLabels = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::CentFT0Cs, aod::CentNTPVs>;
92105

93106
PresliceUnsorted<CollisionWLabels> colPerMcCollision = aod::mccollisionlabel::mcCollisionId;
94107
Preslice<TrackWLabels> tracksPerCollision = aod::track::collisionId;
108+
Preslice<aod::McParticles> perMcCollision = aod::mcparticle::mcCollisionId;
95109

96110
HistogramRegistry registry{"registry", {}};
97-
std::shared_ptr<TH1> hCharmPerCollImpPar;
111+
std::shared_ptr<TH1> hCharmPerCollImpPar, hCollisions;
112+
113+
AxisSpec impParBins = {200, 0, 20};
114+
AxisSpec deltaXYbins = {200, -0.05, 0.05};
115+
AxisSpec deltaZbins = {200, -10, 10};
98116

99117
bool isCharm(int pdg)
100118
{
@@ -105,7 +123,7 @@ struct checkMcPvContr {
105123
return false;
106124
}
107125

108-
bool isBeauty(int pdg)
126+
bool isBeauty(int pdg) // if needed in the future
109127
{
110128
if (std::abs(pdg) / 1000 == 5)
111129
return true;
@@ -116,49 +134,92 @@ struct checkMcPvContr {
116134

117135
void init(InitContext&)
118136
{
119-
registry.add("hCharmImpPar", ";Impact parameter (fm);Charm counts", {HistType::kTH1F, {{200, 0, 20}}});
120-
registry.add("hCollImpPar", ";Impact parameter (fm);Counts", {HistType::kTH1F, {{200, 0, 20}}});
121-
hCharmPerCollImpPar = registry.add<TH1>("hCharmPerCollImpPar", ";Impact parameter (fm);Charm counts per collision", {HistType::kTH1F, {{200, 0, 20}}});
137+
registry.add("hCharmImpPar", ";Impact parameter (fm);Charm counts", {HistType::kTH1F, {impParBins}});
138+
registry.add("hCollImpPar", ";Impact parameter (fm);Counts", {HistType::kTH1F, {impParBins}});
139+
hCharmPerCollImpPar = registry.add<TH1>("hCharmPerCollImpPar", ";Impact parameter (fm);Charm counts per collision", {HistType::kTH1F, {impParBins}});
122140

123-
registry.add("hDeltaX", ";#DeltaX (cm);Counts", {HistType::kTH1F, {{200, -0.01, 0.01}}});
124-
registry.add("hDeltaY", ";#DeltaY (cm);Counts", {HistType::kTH1F, {{200, -0.01, 0.01}}});
125-
registry.add("hDeltaZ", ";#DeltaZ (cm);Counts", {HistType::kTH1F, {{200, -0.01, 0.01}}});
141+
registry.add("hDeltaX", ";#DeltaX (cm);Counts", {HistType::kTH1F, {{deltaXYbins}}});
142+
registry.add("hDeltaY", ";#DeltaY (cm);Counts", {HistType::kTH1F, {{deltaXYbins}}});
143+
registry.add("hDeltaZ", ";#DeltaZ (cm);Counts", {HistType::kTH1F, {{deltaZbins}}});
144+
145+
registry.add("hDeltaX_NPV_lt2000", ";#DeltaX (cm);Counts", {HistType::kTH1F, {{deltaXYbins}}});
146+
registry.add("hDeltaY_NPV_lt2000", ";#DeltaY (cm);Counts", {HistType::kTH1F, {{deltaXYbins}}});
147+
registry.add("hDeltaZ_NPV_lt2000", ";#DeltaZ (cm);Counts", {HistType::kTH1F, {{deltaZbins}}});
148+
149+
registry.add("hDeltaX_NPV_gt2000", ";#DeltaX (cm);Counts", {HistType::kTH1F, {{deltaXYbins}}});
150+
registry.add("hDeltaY_NPV_gt2000", ";#DeltaY (cm);Counts", {HistType::kTH1F, {{deltaXYbins}}});
151+
registry.add("hDeltaZ_NPV_gt2000", ";#DeltaZ (cm);Counts", {HistType::kTH1F, {{deltaZbins}}});
152+
153+
registry.add("hDeltaXSngBkg", ";#DeltaX (signal/bkg) (cm);Counts", {HistType::kTH1F, {{200, -10, 10}}});
154+
registry.add("hDeltaYSngBkg", ";#DeltaY (signal/bkg) (cm);Counts", {HistType::kTH1F, {{200, -10, 10}}});
155+
registry.add("hDeltaZSngBkg", ";#DeltaZ (signal/bkg) (cm);Counts", {HistType::kTH1F, {{200, -20, 20}}});
156+
157+
hCollisions = registry.add<TH1>("hCollisions", ";;Counts", {HistType::kTH1F, {{2, 0.5, 2.5}}});
158+
hCollisions->GetXaxis()->SetBinLabel(1, "Generated");
159+
hCollisions->GetXaxis()->SetBinLabel(2, "Reconstructed");
126160
}
127161

128162
void process(CollisionWLabels const& collisions,
129163
TrackWLabels const& tracks,
130164
aod::McParticles const& mcParticles,
131165
aod::McCollisions const& mcCollisions)
132166
{
133-
assert(isCharm(413));
134167
int splitColls{0};
135168
for (const auto& mcColl : mcCollisions) {
136-
const auto collSlice = collisions.sliceBy(colPerMcCollision, mcColl.globalIndex());
137-
int64_t idxCollMostPVContrib{-1};
138-
int nPVContribMost{0};
139-
// First we find the collision with most PV contributors
140-
for (const auto& collision : collSlice) {
141-
if (collision.centFT0C() < 20.f) {
142-
if (collision.numContrib() > nPVContribMost) {
143-
nPVContribMost = collision.numContrib();
144-
idxCollMostPVContrib = collision.globalIndex();
145-
}
169+
registry.fill(HIST("hCollImpPar"), mcColl.impactParameter());
170+
const auto mcPartColl = mcParticles.sliceBy(perMcCollision, mcColl.globalIndex());
171+
double xAvgSgn{0.}, yAvgSgn{0.}, zAvgSgn{0.};
172+
double xAvgBkg{0.}, yAvgBkg{0.}, zAvgBkg{0.};
173+
int nSgn{0}, nBkg{0};
174+
for (const auto& mcPart : mcPartColl) {
175+
if (isCharm(mcPart.pdgCode())) { // charm hadron
176+
registry.fill(HIST("hCharmImpPar"), mcColl.impactParameter());
177+
}
178+
if (mcPart.fromBackgroundEvent()) {
179+
xAvgBkg += mcPart.vx();
180+
yAvgBkg += mcPart.vy();
181+
zAvgBkg += mcPart.vz();
182+
nBkg++;
183+
tracksInj(mcPart.mcCollisionId(), mcPart.vx(), mcPart.vy(), mcPart.vz(), false);
184+
} else {
185+
xAvgSgn += mcPart.vx();
186+
yAvgSgn += mcPart.vy();
187+
zAvgSgn += mcPart.vz();
188+
nSgn++;
189+
tracksInj(mcPart.mcCollisionId(), mcPart.vx(), mcPart.vy(), mcPart.vz(), true);
146190
}
147191
}
192+
registry.fill(HIST("hDeltaXSngBkg"), xAvgSgn / nSgn - xAvgBkg / nBkg);
193+
registry.fill(HIST("hDeltaYSngBkg"), yAvgSgn / nSgn - yAvgBkg / nBkg);
194+
registry.fill(HIST("hDeltaZSngBkg"), zAvgSgn / nSgn - zAvgBkg / nBkg);
195+
196+
const auto collSlice = collisions.sliceBy(colPerMcCollision, mcColl.globalIndex());
148197

149198
// Then we fill the histogram with the distances of the collisions
150199
for (const auto& collision : collSlice) {
151200
const auto collTracks = tracks.sliceBy(tracksPerCollision, collision.globalIndex());
152201
std::vector<int> charmIds{};
153202
int fromSignalEv{0};
203+
if (collision.centFT0C() < 20.f) {
204+
registry.fill(HIST("hDeltaX"), collision.posX() - collision.mcCollision().posX());
205+
registry.fill(HIST("hDeltaY"), collision.posY() - collision.mcCollision().posY());
206+
registry.fill(HIST("hDeltaZ"), collision.posZ() - collision.mcCollision().posZ());
207+
208+
if (collision.numContrib() > 2000) {
209+
registry.fill(HIST("hDeltaX_NPV_gt2000"), collision.posX() - collision.mcCollision().posX());
210+
registry.fill(HIST("hDeltaY_NPV_gt2000"), collision.posY() - collision.mcCollision().posY());
211+
registry.fill(HIST("hDeltaZ_NPV_gt2000"), collision.posZ() - collision.mcCollision().posZ());
212+
} else {
213+
registry.fill(HIST("hDeltaX_NPV_lt2000"), collision.posX() - collision.mcCollision().posX());
214+
registry.fill(HIST("hDeltaY_NPV_lt2000"), collision.posY() - collision.mcCollision().posY());
215+
registry.fill(HIST("hDeltaZ_NPV_lt2000"), collision.posZ() - collision.mcCollision().posZ());
216+
}
217+
}
154218
for (const auto& track : collTracks) {
155219
if (track.has_mcParticle()) {
156220
auto mcPart = track.mcParticle_as<aod::McParticles>();
157221
for (const auto& mother : mcPart.mothers_as<aod::McParticles>()) {
158222
if (isCharm(mother.pdgCode())) { // charm hadron
159-
registry.fill(HIST("hDeltaX"), collision.posX() - collisions.rawIteratorAt(idxCollMostPVContrib).posX());
160-
registry.fill(HIST("hDeltaY"), collision.posY() - collisions.rawIteratorAt(idxCollMostPVContrib).posY());
161-
registry.fill(HIST("hDeltaZ"), collision.posZ() - collisions.rawIteratorAt(idxCollMostPVContrib).posZ());
162223
if (std::find(charmIds.begin(), charmIds.end(), mother.globalIndex()) == charmIds.end()) {
163224
charmIds.push_back(mother.globalIndex());
164225
fromSignalEv += static_cast<int>(!mother.fromBackgroundEvent());
@@ -169,33 +230,20 @@ struct checkMcPvContr {
169230
}
170231
}
171232
checkInj(
172-
mcColl.globalIndex(), mcColl.impactParameter(), mcColl.posX(), mcColl.posY(), mcColl.posZ(), mcColl.t(), collision.collisionTime(),
173-
charmIds.size(), fromSignalEv, collision.numContrib(), collision.centFT0C(), collision.posX(), collision.posY(), collision.posZ(),
174-
collision.bcId());
175-
}
176-
}
177-
for (const auto& mcColl : mcCollisions) {
178-
registry.fill(HIST("hCollImpPar"), mcColl.impactParameter());
179-
}
180-
int count_bkg{0}, count_sgn{0};
181-
for (const auto& mcPart : mcParticles) {
182-
if (isCharm(mcPart.pdgCode())) { // charm hadron
183-
if (!mcPart.fromBackgroundEvent()) {
184-
auto mcCollision = mcPart.mcCollision_as<aod::McCollisions>();
185-
registry.fill(HIST("hCharmImpPar"), mcCollision.impactParameter());
186-
count_sgn++;
187-
} else {
188-
count_bkg++;
189-
}
233+
mcColl.globalIndex(), mcColl.impactParameter(),
234+
mcColl.posX(), mcColl.posY(), mcColl.posZ(), mcColl.t(),
235+
collision.posX(), collision.posY(), collision.posZ(), collision.collisionTime(),
236+
charmIds.size(), fromSignalEv, collision.numContrib(), collision.centFT0C(), collision.bcId());
190237
}
191238
}
239+
192240
hCharmPerCollImpPar->Divide(registry.get<TH1>(HIST("hCharmImpPar")).get(), registry.get<TH1>(HIST("hCollImpPar")).get(), 1, 1);
193-
LOG(info) << "Number of bkgev particles: " << count_bkg;
194-
LOG(info) << "Number of sngev particles: " << count_sgn;
241+
hCollisions->Fill(1, mcCollisions.size());
242+
hCollisions->Fill(2, collisions.size());
195243
}
196244
};
197245

198246
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
199247
{
200-
return WorkflowSpec{adaptAnalysisTask<checkMcPvContr>(cfgc)};
248+
return WorkflowSpec{adaptAnalysisTask<taskMcInjection>(cfgc)};
201249
}

0 commit comments

Comments
 (0)