Skip to content

Commit 7657008

Browse files
atrioloatriolo
committed
Fix nonPromptCascade task
1 parent eb933ce commit 7657008

1 file changed

Lines changed: 109 additions & 61 deletions

File tree

‎PWGLF/Tasks/nonPromptCascade.cxx‎

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

12-
#include "Framework/AnalysisTask.h"
13-
#include "Framework/AnalysisDataModel.h"
14-
#include "Framework/ASoA.h"
15-
#include "Framework/ASoAHelpers.h"
16-
#include "Framework/HistogramRegistry.h"
17-
#include "Framework/runDataProcessing.h"
12+
#include "CCDB/BasicCCDBManager.h"
1813
#include "Common/DataModel/Centrality.h"
1914
#include "Common/DataModel/EventSelection.h"
2015
#include "Common/DataModel/PIDResponse.h"
2116
#include "Common/DataModel/TrackSelectionTables.h"
2217
#include "Common/Core/RecoDecay.h"
2318
#include "Common/Core/trackUtilities.h"
24-
#include "ReconstructionDataFormats/DCA.h"
25-
#include "ReconstructionDataFormats/Track.h"
26-
#include "PWGLF/DataModel/LFStrangenessTables.h"
27-
#include "CCDB/BasicCCDBManager.h"
28-
#include "DetectorsBase/Propagator.h"
2919
#include "DataFormatsParameters/GRPMagField.h"
3020
#include "DataFormatsParameters/GRPObject.h"
3121
#include "DataFormatsTPC/BetheBlochAleph.h"
32-
33-
#include "PWGLF/DataModel/LFNonPromptCascadeTables.h"
34-
22+
#include "DCAFitter/DCAFitterN.h"
23+
#include "DetectorsBase/Propagator.h"
24+
#include "Framework/AnalysisTask.h"
25+
#include "Framework/AnalysisDataModel.h"
26+
#include "Framework/ASoA.h"
27+
#include "Framework/ASoAHelpers.h"
28+
#include "Framework/HistogramRegistry.h"
29+
#include "Framework/runDataProcessing.h"
30+
#include "PWGHF/Core/PDG.h"
31+
#include "PWGLF/DataModel/LFStrangenessTables.h"
32+
#include "ReconstructionDataFormats/DCA.h"
33+
#include "ReconstructionDataFormats/Track.h"
34+
// #include "PWGLF/DataModel/LFNonPromptCascadeTables.h"
3535

3636
using namespace o2;
3737
using namespace o2::framework;
@@ -57,20 +57,27 @@ std::shared_ptr<TH1> invMassBCOmega;
5757
std::shared_ptr<TH1> invMassACOmega;
5858
std::shared_ptr<TH1> invMassBCXi;
5959
std::shared_ptr<TH1> invMassACXi;
60+
std::shared_ptr<TH1> invMassBCV0;
61+
std::shared_ptr<TH1> invMassACV0;
6062

6163
} // namespace
6264

6365
struct NonPromptCascadeTask {
6466

65-
Produces<o2::aod::NonPromptCascadeTable> nonPromptCascadeTable;
66-
Produces<o2::aod::NonPromptCascadeTableMC> nonPromptCascadeTableMC;
67+
// Produces<o2::aod::NonPromptCascadeTable> nonPromptCascadeTable;
68+
// Produces<o2::aod::NonPromptCascadeTableMC> nonPromptCascadeTableMC;
6769

6870
using TracksExtData = soa::Join<aod::TracksIU, aod::TracksCovIU, aod::TracksExtra, aod::pidTPCFullKa, aod::pidTPCFullPi, aod::pidTPCFullPr>;
6971
using TracksExtMC = soa::Join<aod::TracksIU, aod::TracksCovIU, aod::TracksExtra, aod::McTrackLabels, aod::pidTPCFullKa, aod::pidTPCFullPi, aod::pidTPCFullPr>;
7072
using CollisionCandidatesRun3 = soa::Join<aod::Collisions, aod::EvSels>::iterator;
7173

7274
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
73-
Configurable<double> bz{"bz", -50., "magnetic field"};
75+
Configurable<bool> propToDCA{"propToDCA", true, "create tracks version propagated to PCA"};
76+
Configurable<bool> useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"};
77+
Configurable<double> maxR{"maxR", 200., "reject PCA's above this radius"};
78+
Configurable<double> maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"};
79+
Configurable<double> minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"};
80+
Configurable<double> minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations if chi2/chi2old > this"};
7481
Configurable<int> cfgMaterialCorrection{"cfgMaterialCorrection", static_cast<int>(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Type of material correction"};
7582
Configurable<std::string> cfgGRPmagPath{"cfgGRPmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
7683
Configurable<std::string> cfgGRPpath{"cfgGRPpath", "GLO/GRP/GRP", "Path of the grp file"};
@@ -80,7 +87,7 @@ struct NonPromptCascadeTask {
8087

8188
Service<o2::ccdb::BasicCCDBManager> ccdb;
8289
int mRunNumber = 0;
83-
float mBz = 0.f;
90+
float bz = 0.f;
8491

8592
HistogramRegistry registry{
8693
"registry",
@@ -120,27 +127,27 @@ struct NonPromptCascadeTask {
120127
{"h_buildermassvspt_Xi", "Mass (from builder) vs p_{T};Mass (GeV/#it{c}^2);p_{T} (GeV/#it{c})", {HistType::kTH2D, {{125, 1.296, 1.346}, {50, 0., 10.}}}},
121128
{"h_massvsmass_Xi", "Mass vs mass;Mass (GeV/#it{c}^{2});Mass (GeV/#it{c}^{2})", {HistType::kTH2D, {{125, 1.296, 1.346}, {125, 1.296, 1.346}}}},
122129
{"h_bachelorsign_Xi", "Bachelor sign;Sign;Counts", {HistType::kTH1D, {{6, -3., 3.}}}},
130+
131+
{"h_massvspt_V0", "Mass vs p_{T};Mass (GeV/#it{c}^2);p_{T} (GeV/#it{c})", {HistType::kTH2D, {{125, 1.090, 1.140}, {50, 0., 10.}}}},
132+
123133
}};
124134

125135
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
126136
{
127-
128-
if (mRunNumber == bc.runNumber()) {
129-
return;
130-
}
131-
LOG(debug) << "Run number: " << mRunNumber << " bc runNumber: " << bc.runNumber();
132-
133-
auto run3grp_timestamp = bc.timestamp();
134-
mRunNumber = bc.runNumber();
135-
136-
if (o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(cfgGRPpath, run3grp_timestamp)) {
137-
o2::base::Propagator::initFieldFromGRP(grpo);
138-
LOG(debug) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << mBz << " kZG";
139-
} else if (o2::parameters::GRPMagField* grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(cfgGRPmagPath, run3grp_timestamp)) {
140-
o2::base::Propagator::initFieldFromGRP(grpmag);
141-
LOG(debug) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << mBz << " kZG";
142-
} else {
143-
LOG(fatal) << "Got nullptr from CCDB for path " << cfgGRPpath << " of object GRPMagField and " << cfgGRPmagPath << " of object GRPObject for timestamp " << run3grp_timestamp;
137+
if (mRunNumber != bc.runNumber()) {
138+
mRunNumber = bc.runNumber();
139+
auto timestamp = bc.timestamp();
140+
141+
if (o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(cfgGRPpath, timestamp)) {
142+
o2::base::Propagator::initFieldFromGRP(grpo);
143+
bz = grpo->getNominalL3Field();
144+
} else if (o2::parameters::GRPMagField* grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(cfgGRPmagPath, timestamp)) {
145+
o2::base::Propagator::initFieldFromGRP(grpmag);
146+
bz = std::lround(5.f * grpmag->getL3Current() / 30000.f);
147+
LOG(debug)<<"bz = "<<bz;
148+
} else {
149+
LOG(fatal) << "Got nullptr from CCDB for path " << cfgGRPmagPath << " of object GRPMagField and " << cfgGRPpath << " of object GRPObject for timestamp " << timestamp;
150+
}
144151
}
145152
return;
146153
}
@@ -180,44 +187,34 @@ struct NonPromptCascadeTask {
180187
invMassACOmega = registry.add<TH1>("h_invariantmass_afterCuts_Omega", "Invariant Mass (GeV/#it{c}^{2})", HistType::kTH1D, {{125, 1.650, 1.700, "Invariant Mass (GeV/#it{c}^{2})"}});
181188
invMassBCXi = registry.add<TH1>("h_invariantmass_beforeCuts_Xi", "Invariant Mass (GeV/#it{c}^{2})", HistType::kTH1D, {{125, 1.296, 1.346, "Invariant Mass (GeV/#it{c}^{2})"}});
182189
invMassACXi = registry.add<TH1>("h_invariantmass_afterCuts_Xi", "Invariant Mass (GeV/#it{c}^{2})", HistType::kTH1D, {{125, 1.296, 1.346, "Invariant Mass (GeV/#it{c}^{2})"}});
183-
}
190+
invMassBCV0 = registry.add<TH1>("h_invariantmass_beforeCuts_V0", "Invariant Mass (GeV/#it{c}^{2})", HistType::kTH1D, {{125, 1.090, 1.140, "Invariant Mass (GeV/#it{c}^{2})"}});
191+
invMassACV0 = registry.add<TH1>("h_invariantmass_afterCuts_V0", "Invariant Mass (GeV/#it{c}^{2})", HistType::kTH1D, {{125, 1.090, 1.140, "Invariant Mass (GeV/#it{c}^{2})"}});
192+
}
184193

185194
template <typename TC, typename T, typename B, typename PR, typename PI>
186195
void fillCascadeDCA(TC const& trackedCascade, T const track, B const bachelor, PR const& protonTrack, PI const& pionTrack, o2::dataformats::VertexBase primaryVertex, bool isOmega)
187196
{
188197
const auto matCorr = static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value);
189-
std::array<std::array<float, 3>, 2> momenta;
190-
std::array<double, 2> masses;
191198
auto trackCovTrk = getTrackParCov(track);
192199
o2::dataformats::DCA impactParameterTrk;
193200

194-
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovTrk, mBz, 2.f, matCorr, &impactParameterTrk)) {
195-
momenta[0] = {protonTrack.px() + pionTrack.px(), protonTrack.py() + pionTrack.py(), protonTrack.pz() + pionTrack.pz()};
196-
momenta[1] = {bachelor.px(), bachelor.py(), bachelor.pz()};
197-
masses = {constants::physics::MassLambda, constants::physics::MassKaonCharged};
198-
const auto massOmega = RecoDecay::m(momenta, masses);
199-
201+
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovTrk, bz, 2.f, matCorr, &impactParameterTrk)) {
200202
if (protonTrack.hasTPC() && pionTrack.hasTPC()) {
201203
if (isOmega) {
202204
registry.fill(HIST("h_dca_Omega"), TMath::Sqrt(impactParameterTrk.getR2()));
203205
registry.fill(HIST("h_dcaxy_Omega"), impactParameterTrk.getY());
204206
registry.fill(HIST("h_dcaz_Omega"), impactParameterTrk.getZ());
205207
registry.fill(HIST("h_dcavspt_Omega"), impactParameterTrk.getY(), track.pt());
206208
registry.fill(HIST("h_dcavsr_Omega"), impactParameterTrk.getY(), std::hypot(track.x(), track.y()));
207-
registry.fill(HIST("h_massvspt_Omega"), massOmega, track.pt());
208209
}
209210
}
210211

211-
masses = {constants::physics::MassLambda, constants::physics::MassPionCharged};
212-
const auto massXi = RecoDecay::m(momenta, masses);
213-
214212
if (protonTrack.hasTPC() && pionTrack.hasTPC()) {
215213
registry.fill(HIST("h_dca_Xi"), TMath::Sqrt(impactParameterTrk.getR2()));
216214
registry.fill(HIST("h_dcaxy_Xi"), impactParameterTrk.getY());
217215
registry.fill(HIST("h_dcaz_Xi"), impactParameterTrk.getZ());
218216
registry.fill(HIST("h_dcavspt_Xi"), impactParameterTrk.getY(), track.pt());
219217
registry.fill(HIST("h_dcavsr_Xi"), impactParameterTrk.getY(), std::hypot(track.x(), track.y()));
220-
registry.fill(HIST("h_massvspt_Xi"), massXi, track.pt());
221218
}
222219
}
223220
}
@@ -230,7 +227,7 @@ struct NonPromptCascadeTask {
230227

231228
auto trackCovBach = getTrackParCov(bachelor);
232229
o2::dataformats::DCA impactParameterBach;
233-
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovBach, mBz, 2.f, matCorr, &impactParameterBach)) {
230+
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovBach, bz, 2.f, matCorr, &impactParameterBach)) {
234231
if (isOmega) {
235232
if (bachelor.sign() < 0) {
236233
registry.fill(HIST("h_bachdcaxyM_Omega"), impactParameterBach.getY());
@@ -257,7 +254,7 @@ struct NonPromptCascadeTask {
257254

258255
auto trackCovNtrack = getTrackParCov(pionTrack);
259256
o2::dataformats::DCA impactParameterNtrack;
260-
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovNtrack, mBz, 2.f, matCorr, &impactParameterNtrack)) {
257+
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovNtrack, bz, 2.f, matCorr, &impactParameterNtrack)) {
261258
if (isOmega) {
262259
registry.fill(HIST("h_ntrackdcavspt_Omega"), impactParameterNtrack.getY(), pionTrack.pt());
263260
}
@@ -266,7 +263,7 @@ struct NonPromptCascadeTask {
266263

267264
auto trackCovPtrack = getTrackParCov(protonTrack);
268265
o2::dataformats::DCA impactParameterPtrack;
269-
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovPtrack, mBz, 2.f, matCorr, &impactParameterPtrack)) {
266+
if (o2::base::Propagator::Instance()->propagateToDCA(primaryVertex, trackCovPtrack, bz, 2.f, matCorr, &impactParameterPtrack)) {
270267
if (isOmega) {
271268
registry.fill(HIST("h_ptrackdcavspt_Omega"), impactParameterPtrack.getY(), protonTrack.pt());
272269
}
@@ -367,10 +364,17 @@ struct NonPromptCascadeTask {
367364
}
368365

369366
registry.fill(HIST("h_PIDcutsXi"), 5, trackedCascade.xiMass());
370-
registry.fill(HIST("h_PIDcutsOmega"), 5, trackedCascade.omegaMass());
367+
368+
if (isOmega){
369+
registry.fill(HIST("h_PIDcutsOmega"), 5, trackedCascade.omegaMass());
370+
invMassACOmega->Fill(trackedCascade.omegaMass());
371+
registry.fill(HIST("h_massvspt_Omega"), trackedCascade.omegaMass(), track.pt());
372+
}
373+
374+
registry.fill(HIST("h_PIDcutsXi"), 5, trackedCascade.xiMass());
371375

372376
invMassACXi->Fill(trackedCascade.xiMass());
373-
invMassACOmega->Fill(trackedCascade.omegaMass());
377+
registry.fill(HIST("h_massvspt_Xi"), trackedCascade.xiMass(), track.pt());
374378

375379
fillCascadeDCA(trackedCascade, track, bachelor, protonTrack, pionTrack, primaryVertex, isOmega);
376380

@@ -421,6 +425,15 @@ struct NonPromptCascadeTask {
421425

422426
const auto primaryVertex = getPrimaryVertex(collision);
423427

428+
o2::vertexing::DCAFitterN<2> df2;
429+
df2.setBz(bz);
430+
df2.setPropagateToPCA(propToDCA);
431+
df2.setMaxR(maxR);
432+
df2.setMaxDZIni(maxDZIni);
433+
df2.setMinParamChange(minParamChange);
434+
df2.setMinRelChi2Change(minRelChi2Change);
435+
df2.setUseAbsDCA(useAbsDCA);
436+
424437
for (const auto& trackedCascade : trackedCascades) {
425438

426439
isOmega = false;
@@ -436,20 +449,47 @@ struct NonPromptCascadeTask {
436449

437450
std::array<std::array<float, 3>, 2> momenta;
438451
std::array<double, 2> masses;
439-
momenta[0] = {protonTrack.px() + pionTrack.px(), protonTrack.py() + pionTrack.py(), protonTrack.pz() + pionTrack.pz()};
440-
momenta[1] = {bachelor.px(), bachelor.py(), bachelor.pz()};
441-
masses = {constants::physics::MassLambda, constants::physics::MassKaonCharged};
452+
453+
//track propagation
454+
o2::track::TrackParCov trackParCovV0;
455+
o2::track::TrackPar trackParV0;
456+
o2::track::TrackPar trackParBachelor;
457+
if (df2.process(getTrackParCov(pionTrack), getTrackParCov(protonTrack))) {
458+
trackParCovV0 = df2.createParentTrackParCov(0);
459+
if (df2.process(trackParCovV0, getTrackParCov(bachelor))) {
460+
trackParV0 = df2.getTrackParamAtPCA(0);
461+
trackParBachelor = df2.getTrackParamAtPCA(1);
462+
trackParV0.getPxPyPzGlo(momenta[0]); // getting the V0 momentum
463+
trackParBachelor.getPxPyPzGlo(momenta[1]); // getting the bachelor momentum
464+
} else {
465+
continue;
466+
}
467+
} else {
468+
continue;
469+
}
470+
471+
// Omega
472+
masses = {o2::analysis::pdg::MassLambda0, o2::analysis::pdg::MassKPlus};
442473
const auto massOmega = RecoDecay::m(momenta, masses);
443-
masses = {constants::physics::MassLambda, constants::physics::MassPionCharged};
474+
475+
// Xi
476+
masses = {o2::analysis::pdg::MassLambda0, o2::analysis::pdg::MassPiPlus};
444477
const auto massXi = RecoDecay::m(momenta, masses);
445478

479+
// Lambda
480+
masses = {o2::analysis::pdg::MassProton, o2::analysis::pdg::MassPiMinus};
481+
momenta[0] = {protonTrack.px(), protonTrack.py(), protonTrack.pz()};
482+
momenta[1] = {pionTrack.px(), pionTrack.py(), pionTrack.pz()};
483+
const auto v0mass = RecoDecay::m(momenta, masses);
484+
446485
////Omega hypohesis -> rejecting Xi
447486
if (TMath::Abs(massXi - constants::physics::MassXiMinus) > 0.005) {
448487
isOmega = true;
449488
invMassBCOmega->Fill(massOmega);
450489
}
451490

452491
invMassBCXi->Fill(massXi);
492+
invMassBCV0->Fill(v0mass);
453493

454494
registry.fill(HIST("h_PIDcutsXi"), 0, massXi);
455495
registry.fill(HIST("h_PIDcutsOmega"), 0, massOmega);
@@ -503,11 +543,19 @@ struct NonPromptCascadeTask {
503543
continue;
504544
}
505545

546+
if (isOmega){
547+
registry.fill(HIST("h_PIDcutsOmega"), 5, massOmega);
548+
invMassACOmega->Fill(massOmega);
549+
registry.fill(HIST("h_massvspt_Omega"), massOmega, track.pt());
550+
}
551+
506552
registry.fill(HIST("h_PIDcutsXi"), 5, massXi);
507-
registry.fill(HIST("h_PIDcutsOmega"), 5, massOmega);
508553

509554
invMassACXi->Fill(massXi);
510-
invMassACOmega->Fill(massOmega);
555+
registry.fill(HIST("h_massvspt_Xi"), massXi, track.pt());
556+
557+
invMassACV0->Fill(v0mass);
558+
registry.fill(HIST("h_massvspt_V0"), v0mass, track.pt());
511559

512560
fillCascadeDCA(trackedCascade, track, bachelor, protonTrack, pionTrack, primaryVertex, isOmega);
513561
fillDauDCA(trackedCascade, bachelor, protonTrack, pionTrack, primaryVertex, isOmega);
@@ -520,4 +568,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
520568
{
521569
return WorkflowSpec{
522570
adaptAnalysisTask<NonPromptCascadeTask>(cfgc)};
523-
}
571+
}

0 commit comments

Comments
 (0)