Skip to content

Commit cc67213

Browse files
lubynetslubynets
committed
add AP V0<-Casc; implement selection of cascades
1 parent 95ab00e commit cc67213

File tree

2 files changed

+204
-49
lines changed

2 files changed

+204
-49
lines changed

‎PWGDQ/Tasks/v0selector.cxx‎

Lines changed: 185 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -103,6 +103,7 @@ struct v0selector {
103103
Configurable<float> cutMassOmegaLow{"cutMassOmegaLow", 1.667, "cutMassOmegaLow"};
104104

105105
Configurable<bool> produceV0ID{"produceV0ID", false, "Produce additional V0ID table"};
106+
Configurable<bool> selectCascades{"selectCascades", false, "Select cascades in addition to v0s"};
106107
Configurable<bool> produceCascID{"produceCascID", false, "Produce additional CascID table"};
107108

108109
enum { // Reconstructed V0
@@ -205,10 +206,21 @@ struct v0selector {
205206
Configurable<int> mincrossedrows{"mincrossedrows", 70, "min crossed rows"};
206207
Configurable<float> maxchi2tpc{"maxchi2tpc", 4.0, "max chi2/NclsTPC"};
207208
Configurable<bool> fillhisto{"fillhisto", false, "flag to fill histograms"};
208-
// cutNsigmaElTPC, cutNsigmaPiTPC, cutNsigmaPrTPC
209+
// Cascade-related Configurables
210+
Configurable<float> cascDcaMax{"cascDcaMax", 0.3, "DCA cascade daughters"};
211+
Configurable<float> cascV0DcaMax{"cascV0DcaMax", 0.3, "DCA V0 daughters of the cascade"};
212+
Configurable<float> cascRadiusMin{"cascRadiusMin", 0.0, "Cascade Radius min"};
213+
Configurable<float> cascRadiusMax{"cascRadiusMax", 90.0, "Cascade Radius max"};
214+
Configurable<float> cascV0RadiusMin{"cascV0RadiusMin", 0.0, "V0 of the Cascade Radius min"};
215+
Configurable<float> cascV0RadiusMax{"cascV0RadiusMax", 90.0, "V0 of the Cascade Radius max"};
216+
Configurable<float> cascCosinePAMin{"cascCosinePAMin", 0.998, "Cascade CosPA min"};
217+
Configurable<float> cascV0CosinePAMin{"cascV0CosinePAMin", 0.995, "V0 of the Cascade CosPA min"};
218+
Configurable<float> cascV0CosinePAMax{"cascV0CosinePAMax", 1.000, "V0 of the Cascade CosPA max"};
219+
// cutNsigmaElTPC, cutNsigmaPiTPC, cutNsigmaPrTPC, cutNsigmaKaTPC
209220
Configurable<float> cutNsigmaElTPC{"cutNsigmaElTPC", 5.0, "cutNsigmaElTPC"};
210221
Configurable<float> cutNsigmaPiTPC{"cutNsigmaPiTPC", 5.0, "cutNsigmaPiTPC"};
211222
Configurable<float> cutNsigmaPrTPC{"cutNsigmaPrTPC", 5.0, "cutNsigmaPrTPC"};
223+
Configurable<float> cutNsigmaKaTPC{"cutNsigmaKaTPC", 5.0, "cutNsigmaKaTPC"};
212224

213225
HistogramRegistry registry{"registry"};
214226
void init(o2::framework::InitContext&)
@@ -236,10 +248,35 @@ struct v0selector {
236248
registry.add("hV0APplot", "hV0APplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}});
237249
registry.add("hV0APplotSelected", "hV0APplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}});
238250
registry.add("hV0Psi", "hV0Psi", HistType::kTH2F, {{100, 0, TMath::PiOver2()}, {100, 0, 0.1}});
239-
registry.add("hCascAPplot", "hCascAPplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {300, 0.0f, 0.3f}});
240-
registry.add("hCascAPplotSelected", "hCascAPplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {300, 0.0f, 0.3f}});
241-
registry.add("hMassOmega", "hMassOmega", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 1.62f, 1.72f}});
242-
registry.add("hMassAntiOmega", "hMassAntiOmega", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 1.62f, 1.72f}});
251+
if (selectCascades) {
252+
registry.add("hCascPt", "pT", HistType::kTH1F, {{100, 0.0f, 10}});
253+
registry.add("hCascEtaPhi", "#eta vs. #varphi", HistType::kTH2F, {{63, 0, 6.3}, {20, -1.0f, 1.0f}});
254+
registry.add("hCascDCAxyPosToPV", "hCascDCAxyPosToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}});
255+
registry.add("hCascDCAxyNegToPV", "hCascDCAxyNegToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}});
256+
registry.add("hCascDCAxyBachToPV", "hCascDCAxyBachToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}});
257+
registry.add("hCascDCAzPosToPV", "hCascDCAzPosToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}});
258+
registry.add("hCascDCAzNegToPV", "hCascDCAzNegToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}});
259+
registry.add("hCascDCAzBachToPV", "hCascDCAzBachToPV", HistType::kTH1F, {{1000, -5.0f, 5.0f}});
260+
registry.add("hCascAPplot", "hCascAPplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {300, 0.0f, 0.3f}});
261+
registry.add("hCascV0APplot", "hCascV0APplot", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}});
262+
registry.add("hCascAPplotSelected", "hCascAPplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {300, 0.0f, 0.3f}});
263+
registry.add("hCascV0APplotSelected", "hCascV0APplotSelected", HistType::kTH2F, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}});
264+
registry.add("hCascRadius", "hCascRadius", HistType::kTH1F, {{1000, 0.0f, 100.0f}});
265+
registry.add("hCascV0Radius", "hCascV0Radius", HistType::kTH1F, {{1000, 0.0f, 100.0f}});
266+
registry.add("hCascCosPA", "hCascCosPA", HistType::kTH1F, {{50, 0.95f, 1.0f}});
267+
registry.add("hCascV0CosPA", "hCascV0CosPA", HistType::kTH1F, {{50, 0.95f, 1.0f}});
268+
registry.add("hMassOmega", "hMassOmega", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 1.62f, 1.72f}});
269+
registry.add("hMassAntiOmega", "hMassAntiOmega", HistType::kTH2F, {{900, 0.0f, 90.0f}, {100, 1.62f, 1.72f}});
270+
registry.add("hCascDCADau", "hCascDCADau", HistType::kTH1F, {{1000, 0.0f, 10.0f}});
271+
registry.add("hCascV0DCADau", "hCascV0DCADau", HistType::kTH1F, {{1000, 0.0f, 10.0f}});
272+
}
273+
}
274+
275+
if (selectCascades == false && produceCascID == true) {
276+
LOGP(error, "produceCascID is available only when selectCascades is enabled");
277+
}
278+
if (cutAPOmegaUp1 < cutAPOmegaDown1) {
279+
LOGP(error, "cutAPOmegaUp1 must be greater than cutAPOmegaDown1");
243280
}
244281
}
245282

@@ -426,64 +463,165 @@ struct v0selector {
426463
}
427464
}
428465

429-
for (const auto& Casc : Cascs) {
430-
if (fillhisto) {
431-
registry.fill(HIST("hCascCandidate"), 1);
432-
}
433-
434-
// topological selections
466+
if (selectCascades) {
467+
for (const auto& casc : Cascs) {
468+
if (fillhisto) {
469+
registry.fill(HIST("hCascCandidate"), 1);
470+
}
435471

436-
if (fillhisto) {
437-
registry.fill(HIST("hCascCandidate"), 2);
438-
}
472+
if (std::fabs(casc.posTrack_as<FullTracksExt>().eta()) > 0.9) {
473+
continue;
474+
}
475+
if (std::fabs(casc.negTrack_as<FullTracksExt>().eta()) > 0.9) {
476+
continue;
477+
}
478+
if (std::fabs(casc.negTrack_as<FullTracksExt>().eta()) > 0.9) {
479+
continue;
480+
}
439481

440-
const float Cascradius = Casc.cascradius();
482+
if (casc.posTrack_as<FullTracksExt>().tpcNClsCrossedRows() < mincrossedrows) {
483+
continue;
484+
}
485+
if (casc.negTrack_as<FullTracksExt>().tpcNClsCrossedRows() < mincrossedrows) {
486+
continue;
487+
}
488+
if (casc.bachelor_as<FullTracksExt>().tpcNClsCrossedRows() < mincrossedrows) {
489+
continue;
490+
}
491+
if (casc.posTrack_as<FullTracksExt>().tpcChi2NCl() > maxchi2tpc) {
492+
continue;
493+
}
494+
if (casc.negTrack_as<FullTracksExt>().tpcChi2NCl() > maxchi2tpc) {
495+
continue;
496+
}
497+
if (casc.bachelor_as<FullTracksExt>().tpcChi2NCl() > maxchi2tpc) {
498+
continue;
499+
}
500+
if (std::fabs(casc.posTrack_as<FullTracksExt>().dcaXY()) < dcamin) {
501+
continue;
502+
}
503+
if (std::fabs(casc.negTrack_as<FullTracksExt>().dcaXY()) < dcamin) {
504+
continue;
505+
}
506+
if (std::fabs(casc.bachelor_as<FullTracksExt>().dcaXY()) < dcamin) {
507+
continue;
508+
}
509+
if (std::fabs(casc.posTrack_as<FullTracksExt>().dcaXY()) > dcamax) {
510+
continue;
511+
}
512+
if (std::fabs(casc.negTrack_as<FullTracksExt>().dcaXY()) > dcamax) {
513+
continue;
514+
}
515+
if (std::fabs(casc.bachelor_as<FullTracksExt>().dcaXY()) > dcamax) {
516+
continue;
517+
}
441518

442-
const float mOmega = Casc.mOmega();
519+
if (casc.posTrack_as<FullTracksExt>().sign() * casc.negTrack_as<FullTracksExt>().sign() > 0) { // reject same sign pair
520+
continue;
521+
}
443522

444-
const float alpha = Casc.alpha();
445-
const float qt = Casc.qtarm();
523+
if (fillhisto) {
524+
registry.fill(HIST("hCascCandidate"), 2);
525+
}
446526

447-
if (fillhisto) {
448-
registry.fill(HIST("hCascAPplot"), alpha, qt);
449-
}
527+
auto collision = casc.collision_as<aod::Collisions>();
528+
const float collisionX = collision.posX();
529+
const float collisionY = collision.posY();
530+
const float collisionZ = collision.posZ();
450531

451-
const int cascid = checkCascade(alpha, qt);
452-
if (cascid < 0) {
453-
continue;
454-
}
455-
if (fillhisto) {
456-
registry.fill(HIST("hCascAPplotSelected"), alpha, qt);
457-
}
532+
const float cascDca = casc.dcacascdaughters(); // NOTE the name of getter is misleading. In case of no-KF this is sqrt(Chi2)
533+
const float cascV0Dca = casc.dcaV0daughters(); // NOTE the name of getter is misleading. In case of kfDoDCAFitterPreMinimV0 this is sqrt(Chi2)
534+
const float cascRadius = casc.cascradius();
535+
const float cascV0Radius = casc.dcaV0daughters();
536+
const float cascCosinePA = casc.casccosPA(collisionX, collisionY, collisionZ);
537+
const float cascV0CosinePA = casc.v0cosPA(collisionX, collisionY, collisionZ);
458538

459-
auto storeCascAddID = [&](auto gix, auto id) {
460-
if (produceCascID.value) {
461-
cascpidmap[gix] = id;
539+
if (cascDca > cascDcaMax) {
540+
continue;
462541
}
463-
};
542+
if (cascV0Dca > cascV0DcaMax) {
543+
continue;
544+
}
545+
if (cascRadius < cascRadiusMin || cascRadius > cascRadiusMax) {
546+
continue;
547+
}
548+
if (cascV0Radius < cascV0RadiusMin || cascV0Radius > cascV0RadiusMax) {
549+
continue;
550+
}
551+
if (cascCosinePA < cascCosinePAMin) {
552+
continue;
553+
}
554+
if (cascV0CosinePA < cascV0CosinePAMin || cascV0CosinePA > cascV0CosinePAMax) {
555+
continue;
556+
}
557+
558+
const float mOmega = casc.mOmega();
559+
const float alpha = casc.alpha();
560+
const float qt = casc.qtarm();
561+
const float v0Alpha = casc.v0Alpha();
562+
const float v0Qt = casc.v0Qtarm();
464563

465-
if (cascid == kOmega) {
466564
if (fillhisto) {
467-
registry.fill(HIST("hMassOmega"), Cascradius, mOmega);
565+
registry.fill(HIST("hCascPt"), casc.pt());
566+
registry.fill(HIST("hCascEtaPhi"), casc.phi(), casc.eta());
567+
registry.fill(HIST("hCascDCAxyPosToPV"), casc.posTrack_as<FullTracksExt>().dcaXY());
568+
registry.fill(HIST("hCascDCAxyNegToPV"), casc.negTrack_as<FullTracksExt>().dcaXY());
569+
registry.fill(HIST("hCascDCAxyBachToPV"), casc.bachelor_as<FullTracksExt>().dcaXY());
570+
registry.fill(HIST("hCascDCAzPosToPV"), casc.posTrack_as<FullTracksExt>().dcaZ());
571+
registry.fill(HIST("hCascDCAzNegToPV"), casc.negTrack_as<FullTracksExt>().dcaZ());
572+
registry.fill(HIST("hCascDCAzBachToPV"), casc.bachelor_as<FullTracksExt>().dcaZ());
573+
registry.fill(HIST("hCascAPplot"), alpha, qt);
574+
registry.fill(HIST("hCascV0APplot"), v0Alpha, v0Qt);
575+
registry.fill(HIST("hCascRadius"), cascRadius);
576+
registry.fill(HIST("hCascV0Radius"), cascV0Radius);
577+
registry.fill(HIST("hCascCosPA"), cascCosinePA);
578+
registry.fill(HIST("hCascV0CosPA"), cascV0CosinePA);
579+
registry.fill(HIST("hCascDCADau"), cascDca);
580+
registry.fill(HIST("hCascV0DCADau"), cascV0Dca);
468581
}
469-
if (cutMassOmegaLow < mOmega && mOmega < cutMassOmegaHigh) {
470-
pidmap[Casc.bachelorId()] |= (uint8_t(1) << kOmega);
471-
storeCascAddID(Casc.globalIndex(), kOmega);
582+
583+
const int cascid = checkCascade(alpha, qt);
584+
const int v0id = checkV0(v0Alpha, v0Qt);
585+
if (cascid < 0) {
586+
continue;
587+
}
588+
if (v0id != kLambda && v0id != kAntiLambda) {
589+
continue;
472590
}
473-
} else if (cascid == kAntiOmega) {
474591
if (fillhisto) {
475-
registry.fill(HIST("hMassAntiOmega"), Cascradius, mOmega);
592+
registry.fill(HIST("hCascAPplotSelected"), alpha, qt);
593+
registry.fill(HIST("hCascV0APplotSelected"), v0Alpha, v0Qt);
594+
}
595+
596+
auto storeCascAddID = [&](auto gix, auto id) {
597+
if (produceCascID.value) {
598+
cascpidmap[gix] = id;
599+
}
600+
};
601+
602+
if (cascid == kOmega) {
603+
if (fillhisto) {
604+
registry.fill(HIST("hMassOmega"), cascRadius, mOmega);
605+
}
606+
if (cutMassOmegaLow < mOmega && mOmega < cutMassOmegaHigh && std::abs(casc.posTrack_as<FullTracksExt>().tpcNSigmaPr()) < cutNsigmaPrTPC && std::abs(casc.negTrack_as<FullTracksExt>().tpcNSigmaPi()) < cutNsigmaPiTPC && std::abs(casc.bachelor_as<FullTracksExt>().tpcNSigmaKa()) < cutNsigmaKaTPC) {
607+
pidmap[casc.bachelorId()] |= (uint8_t(1) << kOmega);
608+
storeCascAddID(casc.globalIndex(), kOmega);
609+
}
610+
} else if (cascid == kAntiOmega) {
611+
if (fillhisto) {
612+
registry.fill(HIST("hMassAntiOmega"), cascRadius, mOmega);
613+
}
614+
if (cutMassOmegaLow < mOmega && mOmega < cutMassOmegaHigh && std::abs(casc.posTrack_as<FullTracksExt>().tpcNSigmaPi()) < cutNsigmaPiTPC && std::abs(casc.negTrack_as<FullTracksExt>().tpcNSigmaPr()) < cutNsigmaPrTPC && std::abs(casc.bachelor_as<FullTracksExt>().tpcNSigmaKa()) < cutNsigmaKaTPC) {
615+
pidmap[casc.bachelorId()] |= (uint8_t(1) << kAntiOmega);
616+
storeCascAddID(casc.globalIndex(), kAntiOmega);
617+
}
476618
}
477-
if (cutMassOmegaLow < mOmega && mOmega < cutMassOmegaHigh) {
478-
pidmap[Casc.bachelorId()] |= (uint8_t(1) << kAntiOmega);
479-
storeCascAddID(Casc.globalIndex(), kAntiOmega);
619+
} // end of Casc loop
620+
if (produceCascID.value) {
621+
for (auto& casc : Cascs) {
622+
cascmapID(cascpidmap[casc.globalIndex()]);
480623
}
481624
}
482-
} // end of Casc loop
483-
if (produceCascID.value) {
484-
for (auto& Casc : Cascs) {
485-
cascmapID(cascpidmap[Casc.globalIndex()]);
486-
}
487625
}
488626
for (auto& track : tracks) {
489627
// printf("setting pidmap[%lld] = %d\n",track.globalIndex(),pidmap[track.globalIndex()]);

‎PWGLF/DataModel/LFStrangenessTables.h‎

Lines changed: 19 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1464,6 +1464,21 @@ DECLARE_SOA_DYNAMIC_COLUMN(PsiPair, psipair, //! psi pair angle
14641464
}
14651465
return std::asin(clipToPM1(argsin));
14661466
});
1467+
1468+
DECLARE_SOA_DYNAMIC_COLUMN(V0Alpha, v0Alpha, //! Armenteros Alpha
1469+
[](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) {
1470+
// No need to divide by momentum of the v0 (as in the v0data namespace) since the ratio of lQl is evaluated
1471+
const float lQlPos = RecoDecay::dotProd(std::array{pxpos, pypos, pzpos}, std::array{pxneg + pxpos, pyneg + pypos, pzneg + pzpos});
1472+
const float lQlNeg = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxneg + pxpos, pyneg + pypos, pzneg + pzpos});
1473+
return (lQlPos - lQlNeg) / (lQlPos + lQlNeg);
1474+
});
1475+
1476+
DECLARE_SOA_DYNAMIC_COLUMN(V0QtArm, v0Qtarm, //! Armenteros Qt
1477+
[](float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) {
1478+
const float momTot2 = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg);
1479+
const float dp = RecoDecay::dotProd(std::array{pxneg, pyneg, pzneg}, std::array{pxpos + pxneg, pypos + pyneg, pzpos + pzneg});
1480+
return std::sqrt(RecoDecay::p2(pxneg, pyneg, pzneg) - dp * dp / momTot2); // qtarm
1481+
});
14671482
} // namespace cascdata
14681483

14691484
//______________________________________________________
@@ -1542,10 +1557,12 @@ DECLARE_SOA_TABLE(StoredCascCores, "AOD", "CASCCORE", //! core information about
15421557
cascdata::BachelorEta<cascdata::PxBach, cascdata::PyBach, cascdata::PzBach>,
15431558
cascdata::BachelorPhi<cascdata::PxBach, cascdata::PyBach>,
15441559

1545-
// Armenteros-Podolanski
1560+
// Armenteros-Podolanski and psi-pair
15461561
cascdata::Alpha<cascdata::PxPos, cascdata::PyPos, cascdata::PzPos, cascdata::PxNeg, cascdata::PyNeg, cascdata::PzNeg, cascdata::PxBach, cascdata::PyBach, cascdata::PzBach, cascdata::Sign>,
15471562
cascdata::QtArm<cascdata::PxPos, cascdata::PyPos, cascdata::PzPos, cascdata::PxNeg, cascdata::PyNeg, cascdata::PzNeg, cascdata::PxBach, cascdata::PyBach, cascdata::PzBach>,
1548-
cascdata::PsiPair<cascdata::PxPos, cascdata::PyPos, cascdata::PzPos, cascdata::PxNeg, cascdata::PyNeg, cascdata::PzNeg, cascdata::PxBach, cascdata::PyBach, cascdata::PzBach, cascdata::Sign>);
1563+
cascdata::PsiPair<cascdata::PxPos, cascdata::PyPos, cascdata::PzPos, cascdata::PxNeg, cascdata::PyNeg, cascdata::PzNeg, cascdata::PxBach, cascdata::PyBach, cascdata::PzBach, cascdata::Sign>,
1564+
cascdata::V0Alpha<cascdata::PxPos, cascdata::PyPos, cascdata::PzPos, cascdata::PxNeg, cascdata::PyNeg, cascdata::PzNeg>,
1565+
cascdata::V0QtArm<cascdata::PxPos, cascdata::PyPos, cascdata::PzPos, cascdata::PxNeg, cascdata::PyNeg, cascdata::PzNeg>);
15491566

15501567
DECLARE_SOA_TABLE(StoredKFCascCores, "AOD", "KFCASCCORE", //!
15511568
cascdata::Sign, cascdata::MXi, cascdata::MOmega,

0 commit comments

Comments
 (0)