Skip to content

Commit 84b1e64

Browse files
AlexianLAlexianL
committed
modifications, not validated
1 parent 0a1c376 commit 84b1e64

1 file changed

Lines changed: 199 additions & 11 deletions

File tree

‎PWGHF/HFC/Tasks/taskFlow.cxx‎

Lines changed: 199 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -22,13 +22,15 @@
2222
#include "PWGHF/DataModel/CandidateReconstructionTables.h"
2323
#include "PWGHF/DataModel/CandidateSelectionTables.h"
2424
#include "PWGHF/DataModel/TrackIndexSkimmingTables.h"
25-
#include "PWGHF/Utils/utilsPid.h"
2625
#include "PWGMM/Mult/DataModel/bestCollisionTable.h"
2726

2827
#include "Common/CCDB/EventSelectionParams.h"
2928
#include "Common/CCDB/RCTSelectionFlags.h"
3029
#include "Common/Core/RecoDecay.h"
30+
#include "Common/Core/fwdtrackUtilities.h"
31+
#include "Common/DataModel/Centrality.h"
3132
#include "Common/DataModel/EventSelection.h"
33+
#include "Common/DataModel/McCollisionExtra.h"
3234
#include "Common/DataModel/Multiplicity.h"
3335
#include "Common/DataModel/TrackSelectionTables.h"
3436

@@ -38,7 +40,6 @@
3840
#include <DataFormatsParameters/GRPMagField.h>
3941
#include <DetectorsCommonDataFormats/AlignParam.h>
4042
#include <FT0Base/Geometry.h>
41-
#include <FV0Base/Geometry.h>
4243
#include <Framework/ASoAHelpers.h>
4344
#include <Framework/AnalysisDataModel.h>
4445
#include <Framework/AnalysisHelpers.h>
@@ -57,7 +58,7 @@
5758
#include <MathUtils/Utils.h>
5859

5960
#include <THn.h>
60-
#include <TPDGCode.h>
61+
#include <TRandom.h>
6162
#include <TString.h>
6263

6364
#include <sys/types.h>
@@ -136,6 +137,7 @@ enum MftTrackSelectionStep {
136137
Pt,
137138
DCAxy,
138139
DCAz,
140+
Chi2OverNdf,
139141
IsLTF,
140142
IsCA,
141143
NMftTrackSelectionSteps
@@ -298,11 +300,13 @@ struct HfTaskFlow {
298300
Configurable<float> etaMftTrackMin{"etaMftTrackMin", -3.36f, "Minimum value for the eta of MFT tracks when used in cut function"};
299301
Configurable<float> etaMftTrackMaxFilter{"etaMftTrackMaxFilter", -2.0f, "Maximum value for the eta of MFT tracks when used in filter"};
300302
Configurable<float> etaMftTrackMinFilter{"etaMftTrackMinFilter", -3.9f, "Minimum value for the eta of MFT tracks when used in filter"};
303+
Configurable<float> maxChi2OverNdf{"maxChi2OverNdf", 1000.f, "maximum chi2/ndf for MFT tracks"};
301304
Configurable<float> mftMaxDCAxy{"mftMaxDCAxy", 2.0f, "Cut on dcaXY for MFT tracks"};
302305
Configurable<float> mftMaxDCAz{"mftMaxDCAz", 2.0f, "Cut on dcaZ for MFT tracks"};
303306
Configurable<int> nClustersMftTrack{"nClustersMftTrack", 5, "Minimum number of clusters for the reconstruction of MFT tracks"};
304307
Configurable<float> ptMftTrackMax{"ptMftTrackMax", 10.0f, "max value of MFT tracks pT when used in cut function"};
305308
Configurable<float> ptMftTrackMin{"ptMftTrackMin", 0.f, "min value of MFT tracks pT when used in cut function"};
309+
Configurable<bool> useMftChi2OverNdfCut{"useMftChi2OverNdfCut", false, "use mft track chi2/ndf cut"};
306310
Configurable<bool> useMftPtCut{"useMftPtCut", false, "if true, use the Mft pt function cut"};
307311
Configurable<bool> useOnlyCATracks{"useOnlyCATracks", false, "if true, use strictly MFT tracks reconstructed with CA algo."};
308312
Configurable<bool> useOnlyLTFTracks{"useOnlyLTFTracks", false, "if true, use strictly MFT tracks reconstructed with LTF algo."};
@@ -342,15 +346,18 @@ struct HfTaskFlow {
342346
using HfCandidatesSelD0 = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0>>;
343347
using HfCandidatesSelLc = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelLc>>;
344348
using FilteredTracksWDcaSel = soa::Filtered<soa::Join<aod::TracksWDca, aod::TrackSelection, aod::TracksExtra>>;
349+
using FilteredTracksWDcaSelWLabels = soa::Filtered<soa::Join<aod::TracksWDca, aod::TrackSelection, aod::TracksExtra, aod::McTrackLabels>>;
345350

346351
using FilteredMftTracks = soa::Filtered<aod::MFTTracks>;
352+
using FilteredMftTracksWCollsMcLabels = soa::Filtered<soa::Join<aod::MFTTracks, aod::McMFTTrackLabels>>;
347353

348354
// =========================
349355
// using declarations : MC
350356
// =========================
351357

352358
using SmallGroupMcCollisions = soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSel, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>;
353359
using FilteredMcCollisionsWMult = soa::Filtered<soa::Join<aod::McCollisions, aod::MultMCExtras>>;
360+
using FilteredMcCollisionsWMultWCollsExtra = soa::Filtered<soa::Join<aod::McCollisions, aod::McCollsExtra, aod::MultMCExtras>>;
354361
using FilteredMcParticles = soa::Filtered<aod::McParticles>;
355362

356363
// =========================
@@ -498,6 +505,7 @@ struct HfTaskFlow {
498505
labelsMftTracksSelection[MftTrackSelectionStep::Pt] = "MFT tracks after pT selection";
499506
labelsMftTracksSelection[MftTrackSelectionStep::DCAxy] = "MFT tracks after DCAxy selection";
500507
labelsMftTracksSelection[MftTrackSelectionStep::DCAz] = "MFT tracks after DCAz selection";
508+
labelsMftTracksSelection[MftTrackSelectionStep::Chi2OverNdf] = "MFT tracks after Chi2OverNdf selection";
501509
labelsMftTracksSelection[MftTrackSelectionStep::IsLTF] = "Linear Track Finder MFT tracks";
502510
labelsMftTracksSelection[MftTrackSelectionStep::IsCA] = "Cellular Automaton MFT tracks";
503511
registry.get<TH1>(HIST("Data/Mft/hMftTracksSelection"))->SetMinimum(0);
@@ -661,8 +669,7 @@ struct HfTaskFlow {
661669
registry.add("Data/hEfficiencyTrigger", "", {HistType::kTH3D, {{configAxis.axisPtTrigger}, {configAxis.axisEtaTrigger}, {configAxis.axisVertex}}});
662670
registry.add("Data/hEfficiencyAssociated", "", {HistType::kTH3D, {{configAxis.axisPtAssoc}, {configAxis.axisEtaAssociated}, {configAxis.axisVertex}}});
663671

664-
registry.add("Data/hMultiplicity_uncorrected", "", {HistType::kTH1D, {configAxis.axisMultiplicity}});
665-
registry.add("Data/hMultiplicity_corrected", "", {HistType::kTH1D, {configAxis.axisMultiplicity}});
672+
registry.add("Data/hMultiplicity_uncorrected_vs_corrected", "", {HistType::kTH2D, {{configAxis.axisMultiplicity}, {configAxis.axisMultiplicity}}});
666673

667674
if (!configTask.doEtaDependentFlow && !configTask.doVariationContainers) {
668675
registry.add("Trig_hist_TPC_MFT", "", {HistType::kTHnSparseF, {{configAxis.axisSamples, configAxis.axisVertex, configAxis.axisPtTrigger}}});
@@ -918,6 +925,13 @@ struct HfTaskFlow {
918925
}
919926
}
920927

928+
if (doprocessTrackEfficiencies) {
929+
registry.add("MC/hEfficiencyTrigger", "", {HistType::kTH3D, {{configAxis.axisPtTrigger}, {configAxis.axisEtaTrigger}, {configAxis.axisVertex}}});
930+
registry.add("MC/hEfficiencyAssociated", "", {HistType::kTH3D, {{configAxis.axisPtTrigger}, {configAxis.axisEtaAssociated}, {configAxis.axisVertex}}});
931+
registry.add("Data/hEfficiencyTrigger", "", {HistType::kTH3D, {{configAxis.axisPtTrigger}, {configAxis.axisEtaTrigger}, {configAxis.axisVertex}}});
932+
registry.add("Data/hEfficiencyAssociated", "", {HistType::kTH3D, {{configAxis.axisPtAssoc}, {configAxis.axisEtaAssociated}, {configAxis.axisVertex}}});
933+
}
934+
921935
} // End of init() function
922936

923937
// =========================
@@ -1434,6 +1448,17 @@ struct HfTaskFlow {
14341448
registry.fill(HIST("Data/Mft/hMftTracksSelection"), MftTrackSelectionStep::DCAz);
14351449
}
14361450

1451+
if (configMft.useMftChi2OverNdfCut) {
1452+
float ndfMftTrack = std::max(2.0f * mftTrack.nClusters() - 5.0f, 1.0f);
1453+
float mftChi2OverNdf = mftTrack.chi2() / ndfMftTrack;
1454+
if (mftChi2OverNdf > configMft.maxChi2OverNdf) {
1455+
return false;
1456+
}
1457+
}
1458+
if (fillHistograms) {
1459+
registry.fill(HIST("Data/Mft/hMftTracksSelection"), MftTrackSelectionStep::Chi2OverNdf);
1460+
}
1461+
14371462
// cut on the track algorithm of MFT tracks
14381463
if (mftTrack.isCA()) {
14391464
if (fillHistograms) {
@@ -1988,16 +2013,33 @@ struct HfTaskFlow {
19882013

19892014
if (configTask.doEtaDependentFlow) {
19902015
if (fitType == isFT0A) {
1991-
registry.fill(HIST("Trig_hist_TPC_FT0A"), sampleIndex, posZ, track1.eta(), triggerWeight); // think about event weight in near future
2016+
if constexpr (std::is_same_v<FilteredMftTracks, TTracksTrig>) {
2017+
registry.fill(HIST("Trig_hist_MFT_FT0A"), sampleIndex, posZ, track1.eta(), triggerWeight); // think about event weight in near future
2018+
} else {
2019+
registry.fill(HIST("Trig_hist_TPC_FT0A"), sampleIndex, posZ, track1.eta(), triggerWeight); // think about event weight in near future
2020+
}
19922021
} else {
1993-
registry.fill(HIST("Trig_hist_TPC_FT0C"), sampleIndex, posZ, track1.eta(), triggerWeight); // think about event weight in near future
2022+
if constexpr (std::is_same_v<FilteredMftTracks, TTracksTrig>) {
2023+
registry.fill(HIST("Trig_hist_MFT_FT0C"), sampleIndex, posZ, track1.eta(), triggerWeight); // think about event weight in near future
2024+
} else {
2025+
registry.fill(HIST("Trig_hist_TPC_FT0C"), sampleIndex, posZ, track1.eta(), triggerWeight); // think about event weight in near future
2026+
}
19942027
}
19952028
}
19962029
if (configTask.doVariationContainers) {
19972030
if (fitType == isFT0A) {
1998-
registry.fill(HIST("Trig_hist_TPC_FT0A"), sampleIndex, posZ, track1.pt(), triggerWeight); // think about event weight in near future
2031+
2032+
if constexpr (std::is_same_v<FilteredMftTracks, TTracksTrig>) {
2033+
registry.fill(HIST("Trig_hist_MFT_FT0A"), sampleIndex, posZ, track1.pt(), triggerWeight); // think about event weight in near future
2034+
} else {
2035+
registry.fill(HIST("Trig_hist_TPC_FT0A"), sampleIndex, posZ, track1.pt(), triggerWeight); // think about event weight in near future
2036+
}
19992037
} else {
2000-
registry.fill(HIST("Trig_hist_TPC_FT0C"), sampleIndex, posZ, track1.pt(), triggerWeight); // think about event weight in near future
2038+
if constexpr (std::is_same_v<FilteredMftTracks, TTracksTrig>) {
2039+
registry.fill(HIST("Trig_hist_MFT_FT0C"), sampleIndex, posZ, track1.pt(), triggerWeight); // think about event weight in near future
2040+
} else {
2041+
registry.fill(HIST("Trig_hist_TPC_FT0C"), sampleIndex, posZ, track1.pt(), triggerWeight); // think about event weight in near future
2042+
}
20012043
}
20022044
}
20032045

@@ -3019,10 +3061,10 @@ struct HfTaskFlow {
30193061
multiplicity = getMultiplicityEstimator(collision, true);
30203062
}
30213063

3022-
registry.fill(HIST("Data/hMultiplicity_uncorrected"), multiplicity);
3064+
auto uncorrected_multiplicity = multiplicity;
3065+
registry.fill(HIST("Data/hMultiplicity_uncorrected_vs_corrected"), uncorrected_multiplicity, multiplicity);
30233066
if (configCollision.useMultiplicityFromTracksCorrected) {
30243067
multiplicity = getCorrectedMultiplicity(tracks);
3025-
registry.fill(HIST("Data/hMultiplicity_corrected"), multiplicity);
30263068
}
30273069

30283070
if (multiplicity < configCollision.minMultiplicity || multiplicity >= configCollision.maxMultiplicity) {
@@ -4274,6 +4316,152 @@ struct HfTaskFlow {
42744316
}
42754317
PROCESS_SWITCH(HfTaskFlow, processMixedMcGen, "MC Gen : Process mixed-event correlations", false);
42764318

4319+
void processTrackEfficiencies(FilteredMcCollisionsWMultWCollsExtra::iterator const& mcCollision,
4320+
SmallGroupMcCollisions const& reconstructedCollisions,
4321+
FilteredTracksWDcaSelWLabels const& tpcTracks, FilteredMftTracksWCollsMcLabels const& /*mftTracks*/,
4322+
soa::Filtered<soa::Join<aod::BestCollisionsFwd3d, aod::McMFTTrackLabels>> const& reassociated3dMftTracks,
4323+
FilteredMcParticles const& mcParticles)
4324+
{
4325+
if (std::abs(mcCollision.posZ()) >= configCollision.zVertexMax) {
4326+
return;
4327+
}
4328+
4329+
// check whether we have a selected reconstructed collision corresponding to the MC collision
4330+
bool hasReconstructedCollision = false;
4331+
for (const auto& reconstructedCollision : reconstructedCollisions) {
4332+
if (!isAcceptedCollision(reconstructedCollision)) {
4333+
continue;
4334+
}
4335+
if (reconstructedCollision.globalIndex() != mcCollision.bestCollisionIndex()) {
4336+
continue;
4337+
}
4338+
// auto groupedTpcTracks = tpcTracks.sliceBy(perColTracks, reconstructedCollision.globalIndex());
4339+
hasReconstructedCollision = true;
4340+
}
4341+
4342+
// fill histogram for MC Gen TPC particles
4343+
for (const auto& particle : mcParticles) {
4344+
auto pdgParticle = pdg->GetParticle(particle.pdgCode());
4345+
// check MC related properties of the particle
4346+
if (!particle.isPhysicalPrimary() || !particle.producedByGenerator()) {
4347+
continue;
4348+
}
4349+
// check charge of the particle
4350+
if (pdgParticle == nullptr || std::abs(pdgParticle->Charge()) == 0) {
4351+
continue;
4352+
}
4353+
// check kinematics of the particle for TPC
4354+
if (std::abs(particle.eta()) <= configTask.etaMcParticlesTriggerMax &&
4355+
particle.pt() >= configTask.ptMcParticlesTriggerMin &&
4356+
particle.pt() <= configTask.ptMcParticlesTriggerMax) {
4357+
if (hasReconstructedCollision) {
4358+
registry.fill(HIST("MC/hEfficiencyTrigger"), mcCollision.posZ(), particle.eta(), particle.pt());
4359+
}
4360+
}
4361+
// check kinematics of the particle for MFT
4362+
if (particle.eta() <= configTask.etaMcParticlesAssocMax &&
4363+
particle.eta() >= configTask.etaMcParticlesAssocMin &&
4364+
particle.pt() >= configTask.ptMcParticlesAssocMin &&
4365+
particle.pt() <= configTask.ptMcParticlesAssocMax) {
4366+
if (hasReconstructedCollision) {
4367+
registry.fill(HIST("MC/hEfficiencyAssociated"), mcCollision.posZ(), particle.eta(), particle.pt());
4368+
}
4369+
}
4370+
}
4371+
4372+
// // fill histogram for MC Gen TPC particles
4373+
// for (const auto& tpcParticle : mcParticles) {
4374+
// auto pdgTpcParticle = pdg->GetParticle(tpcParticle.pdgCode());
4375+
// // check MC related properties of the particle
4376+
// if (!tpcParticle.isPhysicalPrimary() || !tpcParticle.producedByGenerator()) {
4377+
// continue;
4378+
// }
4379+
// // check charge of the particle
4380+
// if (pdgTpcParticle == nullptr || std::abs(pdgTpcParticle->Charge()) == 0) {
4381+
// continue;
4382+
// }
4383+
// // check kinematics of the particle
4384+
// if (std::abs(tpcParticle.eta()) > configTask.etaMcParticlesTriggerMax ||
4385+
// tpcParticle.pt() < configTask.ptMcParticlesTriggerMin ||
4386+
// tpcParticle.pt() > configTask.ptMcParticlesTriggerMax) {
4387+
// continue;
4388+
// }
4389+
// if (hasReconstructedCollision) {
4390+
// registry.fill(HIST("MC/hEfficiencyTrigger"), mcCollision.posZ(), tpcParticle.eta(), tpcParticle.pt());
4391+
// }
4392+
// }
4393+
4394+
// // fill histogram for MC Gen MFT particles
4395+
// for (const auto& mftParticle : mcParticles) {
4396+
// auto pdgMftParticle = pdg->GetParticle(mftParticle.pdgCode());
4397+
// // check MC related properties of the particle
4398+
// if (!mftParticle.isPhysicalPrimary() || !mftParticle.producedByGenerator()) {
4399+
// continue;
4400+
// }
4401+
// // check charge of the particle
4402+
// if (pdgMftParticle == nullptr || std::abs(pdgMftParticle->Charge()) == 0) {
4403+
// continue;
4404+
// }
4405+
// // check kinematics of the particle
4406+
// if (mftParticle.eta() > configTask.etaMcParticlesAssocMax ||
4407+
// mftParticle.eta() < configTask.etaMcParticlesAssocMin ||
4408+
// mftParticle.pt() < configTask.ptMcParticlesAssocMin ||
4409+
// mftParticle.pt() > configTask.ptMcParticlesAssocMax) {
4410+
// continue;
4411+
// }
4412+
4413+
// if (hasReconstructedCollision) {
4414+
// registry.fill(HIST("MC/hEfficiencyAssociated"), mcCollision.posZ(), mftParticle.eta(), mftParticle.pt());
4415+
// }
4416+
// }
4417+
4418+
// fill histogram for reconstructed TPC tracks
4419+
for (const auto& reconstructedCollision : reconstructedCollisions) {
4420+
if (!isAcceptedCollision(reconstructedCollision)) {
4421+
continue;
4422+
}
4423+
if (reconstructedCollision.globalIndex() != mcCollision.bestCollisionIndex()) {
4424+
continue;
4425+
}
4426+
4427+
auto groupedTpcTracks = tpcTracks.sliceBy(perColTracks, reconstructedCollision.globalIndex());
4428+
for (const auto& tpcTrack : groupedTpcTracks) {
4429+
if (!isAcceptedCentralTrack(tpcTrack)) {
4430+
continue;
4431+
}
4432+
if (!tpcTrack.has_mcParticle()) {
4433+
continue;
4434+
}
4435+
auto tpcParticle = tpcTrack.template mcParticle_as<FilteredMcParticles>();
4436+
if (reconstructedCollision.mcCollisionId() != tpcParticle.mcCollisionId()) {
4437+
continue;
4438+
}
4439+
if (tpcParticle.isPhysicalPrimary()) {
4440+
registry.fill(HIST("Data/hEfficiencyTrigger"), mcCollision.posZ(), tpcParticle.eta(), tpcParticle.pt());
4441+
}
4442+
}
4443+
4444+
auto groupedreassociated3dMftTracks = reassociated3dMftTracks.sliceBy(perColReassociated3dTracks, reconstructedCollision.globalIndex());
4445+
for (const auto& reassociated3dMftTrack : reassociated3dMftTracks) {
4446+
if (!reassociated3dMftTrack.has_mcParticle()) {
4447+
continue;
4448+
}
4449+
auto templatedMftTrack = reassociated3dMftTrack.template mfttrack_as<FilteredMftTracksWCollsMcLabels>();
4450+
if (!isAcceptedMftTrack(templatedMftTrack, reassociated3dMftTrack.bestDCAXY(), reassociated3dMftTrack.bestDCAZ(), false)) {
4451+
continue;
4452+
}
4453+
auto mftParticle = templatedMftTrack.template mcParticle_as<FilteredMcParticles>();
4454+
if (reconstructedCollision.mcCollisionId() != mftParticle.mcCollisionId()) {
4455+
continue;
4456+
}
4457+
if (mftParticle.isPhysicalPrimary()) {
4458+
registry.fill(HIST("Data/hEfficiencyAssociated"), mcCollision.posZ(), mftParticle.eta(), mftParticle.pt());
4459+
}
4460+
}
4461+
}
4462+
}
4463+
PROCESS_SWITCH(HfTaskFlow, processTrackEfficiencies, "process track efficiencies for TPC and MFT", false);
4464+
42774465
}; // End of struct
42784466

42794467
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)