Skip to content

Commit 60d431b

Browse files
[PWGLF] added DCAxy p-dependent cut and added rotational method (#15274)
1 parent 022e2b4 commit 60d431b

File tree

1 file changed

+127
-37
lines changed

1 file changed

+127
-37
lines changed

PWGLF/Tasks/Resonances/lambda1520_PbPb.cxx

Lines changed: 127 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,6 @@
2929
#include <Framework/HistogramRegistry.h>
3030
#include <Framework/HistogramSpec.h>
3131

32-
#include <TLorentzVector.h>
33-
#include <TRandom.h>
34-
3532
#include <fairlogger/Logger.h>
3633

3734
using namespace o2;
@@ -56,7 +53,15 @@ struct lambdaAnalysis_pb {
5653
Configurable<float> cPMin{"cPMin", 0., "Minimum Track p"};
5754
Configurable<float> cEtaCut{"cEtaCut", 0.8, "Pseudorapidity cut"};
5855
Configurable<float> cDcaz{"cDcazMin", 1., "Minimum DCAz"};
59-
Configurable<float> cDcaxy{"cDcaxyMin", 0.1, "Minimum DCAxy"};
56+
57+
Configurable<std::vector<float>> cDcaPtBinsPr{"cDcaPtBinsPr", {0.0f, 0.5f, 1.0f, 2.0f, 3.0f, 5.0f, 1000.0f}, "Proton pT bin edges for DCAxy cut"};
58+
59+
Configurable<std::vector<float>> cDcaXYBinsPr{"cDcaXYBinsPr", {0.020f, 0.015f, 0.010f, 0.007f, 0.005f, 0.004f}, "Proton max |DCAxy| per pT bin (cm)"};
60+
61+
// Kaon DCAxy — pT binned
62+
Configurable<std::vector<float>> cDcaPtBinsKa{"cDcaPtBinsKa", {0.0f, 0.3f, 0.6f, 1.0f, 2.0f, 1000.0f}, "Kaon pT bin edges for DCAxy cut"};
63+
64+
Configurable<std::vector<float>> cDcaXYBinsKa{"cDcaXYBinsKa", {0.025f, 0.018f, 0.012f, 0.008f, 0.004f}, "Kaon max |DCAxy| per pT bin (cm)"};
6065
Configurable<bool> isonlyQC{"isonlyQC", false, "only QC"};
6166
Configurable<bool> isDeepAngle{"isDeepAngle", false, "Deep Angle cut"};
6267
Configurable<double> cfgDeepAngle{"cfgDeepAngle", 0.04, "Deep Angle cut value"};
@@ -108,6 +113,16 @@ struct lambdaAnalysis_pb {
108113
ConfigurableAxis cMixMultBins{"cMixMultBins", {VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 200.0f}, "Mixing bins - multiplicity"};
109114
ConfigurableAxis cMixEPAngle{"cMixEPAngle", {VARIABLE_WIDTH, -1.5708f, -1.25664f, -0.942478f, -0.628319f, 0.f, 0.628319f, 0.942478f, 1.25664f, 1.5708f}, "event plane"};
110115
ConfigurableAxis occupancy_bins{"occupancy_bins", {VARIABLE_WIDTH, 0.0, 100, 500, 600, 1000, 1100, 1500, 1600, 2000, 2100, 2500, 2600, 3000, 3100, 3500, 3600, 4000, 4100, 4500, 4600, 5000, 5100, 9999}, "Binning of the occupancy axis"};
116+
Configurable<int> cNofRotations{"cNofRotations", 10, "Number of rotations for rotational background"};
117+
Configurable<float> rotationalcut{"rotationalcut", 6.f, "Rotational background angle window: PI/rotationalcut"};
118+
119+
// ── MC Event Selection Configurables ─────────────────────────────────────
120+
Configurable<bool> cEvtMCAfterAllCuts{"cEvtMCAfterAllCuts", false, "MC event sel: isInAfterAllCuts"};
121+
Configurable<bool> cEvtMCINELgt0{"cEvtMCINELgt0", false, "MC event sel: isINELgt0"};
122+
Configurable<bool> cEvtMCSel8{"cEvtMCSel8", false, "MC event sel: isInSel8"};
123+
Configurable<bool> cEvtMCVtxIn10{"cEvtMCVtxIn10", false, "MC event sel: isVtxIn10"};
124+
Configurable<bool> cEvtMCTriggerTVX{"cEvtMCTriggerTVX", false, "MC event sel: isTriggerTVX"};
125+
Configurable<bool> cEvtMCRecINELgt0{"cEvtMCRecINELgt0", false, "MC event sel: isRecINELgt0"};
111126
// Histogram Registry.
112127
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
113128

@@ -131,6 +146,7 @@ struct lambdaAnalysis_pb {
131146
AxisSpec axisDCAz = {cDCAzBins, "DCA_{z} (cm)"};
132147

133148
histos.add("Event/h1d_ft0_mult_percentile", "FT0 (%)", kTH2F, {axisCent, axisOccupancy});
149+
histos.add("Event/h_ft0_vz", "Collision Vertex Z position", kTH1F, {{100, -15., 15.}});
134150
if (doprocessMix || doprocessMixDF || doprocessMixepDF) {
135151
histos.add("Event/mixing_vzVsmultpercentile", "FT0(%)", kTH3F, {axisCent, axisVz, axisEP});
136152
}
@@ -208,19 +224,58 @@ struct lambdaAnalysis_pb {
208224
if (std::abs(track.eta()) > cEtaCut)
209225
return false;
210226

211-
if (std::abs(track.dcaZ()) > cDcaz)
227+
if (cPrimaryTrack && !track.isPrimaryTrack())
212228
return false;
213229

214-
if (std::abs(track.dcaXY()) > cDcaxy)
230+
if (cGlobalWoDCATrack && !track.isGlobalTrackWoDCA())
215231
return false;
216232

217-
if (cPrimaryTrack && !track.isPrimaryTrack())
233+
if (cPVContributor && !track.isPVContributor())
218234
return false;
219235

220-
if (cGlobalWoDCATrack && !track.isGlobalTrackWoDCA())
236+
return true;
237+
}
238+
// ── Proton DCA Selection ──────────────────────────────────────────────────
239+
template <typename T>
240+
bool dcaSelectionProton(T const& track, float p)
241+
{
242+
auto ptBinsPr = static_cast<std::vector<float>>(cDcaPtBinsPr);
243+
auto dcaXYPr = static_cast<std::vector<float>>(cDcaXYBinsPr);
244+
int nBinsPr = static_cast<int>(ptBinsPr.size()) - 1;
245+
246+
bool dcaXYPassed = false;
247+
for (int i = 0; i < nBinsPr; i++) {
248+
if (p >= ptBinsPr[i] && p < ptBinsPr[i + 1] &&
249+
std::abs(track.dcaXY()) < dcaXYPr[i])
250+
dcaXYPassed = true;
251+
}
252+
if (!dcaXYPassed)
221253
return false;
222254

223-
if (cPVContributor && !track.isPVContributor())
255+
if (std::abs(track.dcaZ()) > cDcaz)
256+
return false;
257+
258+
return true;
259+
}
260+
261+
// ── Kaon DCA Selection ────────────────────────────────────────────────────
262+
template <typename T>
263+
bool dcaSelectionKaon(T const& track, float p)
264+
{
265+
auto ptBinsKa = static_cast<std::vector<float>>(cDcaPtBinsKa);
266+
auto dcaXYKa = static_cast<std::vector<float>>(cDcaXYBinsKa);
267+
int nBinsKa = static_cast<int>(ptBinsKa.size()) - 1;
268+
269+
bool dcaXYPassed = false;
270+
for (int i = 0; i < nBinsKa; i++) {
271+
if (p >= ptBinsKa[i] && p < ptBinsKa[i + 1] &&
272+
std::abs(track.dcaXY()) < dcaXYKa[i])
273+
dcaXYPassed = true;
274+
}
275+
if (!dcaXYPassed)
276+
return false;
277+
278+
if (std::abs(track.dcaZ()) > cDcaz)
224279
return false;
225280

226281
return true;
@@ -362,8 +417,6 @@ struct lambdaAnalysis_pb {
362417
void fillDataHistos(trackType const& trk1, trackType const& trk2, float mult, int occup = 100)
363418
{
364419

365-
TLorentzVector p1, p2, p;
366-
TRandom* rn = new TRandom();
367420
float p_ptot = 0., k_ptot = 0.;
368421

369422
for (auto const& [trkPr, trkKa] : soa::combinations(soa::CombinationsFullIndexPolicy(trk1, trk2))) {
@@ -412,6 +465,9 @@ struct lambdaAnalysis_pb {
412465
continue;
413466
if (!selectionPIDProton(trkPr, p_ptot) || !selectionPIDKaon(trkKa, k_ptot))
414467
continue;
468+
if (!dcaSelectionProton(trkPr, p_ptot) || !dcaSelectionKaon(trkKa, k_ptot))
469+
continue;
470+
415471
if (isDeepAngle && TMath::ACos((trkPr.pt() * trkKa.pt() + _pzPr * _pzKa) / (p_ptot * k_ptot)) < cfgDeepAngle)
416472
continue;
417473

@@ -461,39 +517,56 @@ struct lambdaAnalysis_pb {
461517

462518
if (isonlyQC)
463519
continue;
464-
// Invariant mass reconstruction.
465-
p1.SetXYZM(_pxPr, _pyPr, _pzPr, MassProton);
466-
p2.SetXYZM(_pxKa, _pyKa, _pzKa, MassKaonCharged);
467-
p = p1 + p2;
468520

469-
if (std::abs(p.Rapidity()) > 0.5)
470-
continue;
521+
std::array<float, 3> pvec0 = {_pxPr, _pyPr, _pzPr};
522+
std::array<float, 3> pvec1 = {_pxKa, _pyKa, _pzKa};
523+
std::array<std::array<float, 3>, 2> arrMomrec = {pvec0, pvec1};
524+
525+
float _M = RecoDecay::m(arrMomrec, std::array{MassProton, MassKaonCharged});
526+
float _pt = RecoDecay::pt(std::array{_pxPr + _pxKa, _pyPr + _pyKa});
471527

472-
auto _M = p.M();
473-
auto _pt = p.Pt();
528+
if (std::abs(RecoDecay::y(std::array{_pxPr + _pxKa, _pyPr + _pyKa, _pzPr + _pzKa}, MassLambda1520)) > 0.5)
529+
continue;
474530

475531
// Apply kinematic cuts.
476-
if (cKinCuts) {
477-
TVector3 v1(_pxPr, _pyPr, _pzPr);
478-
TVector3 v2(_pxKa, _pyKa, _pzKa);
479-
float alpha = v1.Angle(v2);
480-
if (alpha > 1.4 && alpha < 2.4)
481-
continue;
482-
}
483532

484533
// Fill Invariant Mass Histograms.
485534
if constexpr (!mix && !mc) {
486535
if (trkPr.sign() * trkKa.sign() < 0) {
536+
487537
if (trkPr.sign() > 0)
488538
histos.fill(HIST("Analysis/h4d_lstar_invm_US_PM"), _M, _pt, mult, occup);
489539
else
490540
histos.fill(HIST("Analysis/h4d_lstar_invm_US_MP"), _M, _pt, mult, occup);
541+
491542
if (doRotate) {
492-
float theta = rn->Uniform(1.56, 1.58);
493-
p1.RotateZ(theta);
494-
p = p1 + p2;
495-
if (std::abs(p.Rapidity()) < 0.5) {
496-
histos.fill(HIST("Analysis/h4d_lstar_invm_rot"), p.M(), p.Pt(), mult, occup);
543+
for (int i = 0; i < cNofRotations; i++) {
544+
545+
float delta = o2::constants::math::PI / rotationalcut;
546+
float theta2 = (o2::constants::math::PI - delta) + i * (2.f * delta / (cNofRotations - 1));
547+
548+
float phiRot = trkKa.phi() + theta2;
549+
if (phiRot > o2::constants::math::TwoPI)
550+
phiRot -= o2::constants::math::TwoPI;
551+
if (phiRot < 0.f)
552+
phiRot += o2::constants::math::TwoPI;
553+
554+
float _pxKaRot = trkKa.pt() * std::cos(phiRot);
555+
float _pyKaRot = trkKa.pt() * std::sin(phiRot);
556+
557+
std::array<float, 3> pvec0rot = {_pxPr, _pyPr, _pzPr};
558+
std::array<float, 3> pvec1rot = {_pxKaRot, _pyKaRot, _pzKa};
559+
std::array<std::array<float, 3>, 2> arrMomRot = {pvec0rot, pvec1rot};
560+
561+
float _Mrot = RecoDecay::m(arrMomRot, std::array{MassProton, MassKaonCharged});
562+
float _ptRot = RecoDecay::pt(std::array{_pxPr + _pxKaRot, _pyPr + _pyKaRot});
563+
564+
if (std::abs(RecoDecay::y(
565+
std::array{_pxPr + _pxKaRot, _pyPr + _pyKaRot, _pzPr + _pzKa},
566+
MassLambda1520)) > 0.5f)
567+
continue;
568+
569+
histos.fill(HIST("Analysis/h4d_lstar_invm_rot"), _Mrot, _ptRot, mult, occup);
497570
}
498571
}
499572
} else {
@@ -542,24 +615,41 @@ struct lambdaAnalysis_pb {
542615
}
543616

544617
using resoCols = soa::Join<aod::ResoCollisions, aod::ResoEvtPlCollisions>;
618+
using resoMCCols = soa::Join<aod::ResoCollisions, aod::ResoMCCollisions>;
619+
545620
using resoTracks = aod::ResoTracks;
546621

547622
void processData(resoCols::iterator const& collision, resoTracks const& tracks)
548623
{
549624

550625
// LOGF(info, " collisions: Index = %d %d", collision.globalIndex(),tracks.size());
551626
histos.fill(HIST("Event/h1d_ft0_mult_percentile"), collision.cent(), 100);
627+
histos.fill(HIST("Event/h_ft0_vz"), collision.posZ());
628+
552629
fillDataHistos<false, false>(tracks, tracks, collision.cent());
553630
}
554631

555632
PROCESS_SWITCH(lambdaAnalysis_pb, processData, "Process for Same Event Data", true);
556633

557-
void processMC(resoCols::iterator const& collision,
634+
void processMC(resoMCCols::iterator const& collision,
558635
soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& tracks, aod::ResoMCParents const& resoParents)
559636
{
637+
if (cEvtMCAfterAllCuts && !collision.isInAfterAllCuts())
638+
return;
639+
if (cEvtMCINELgt0 && !collision.isINELgt0())
640+
return;
641+
if (cEvtMCSel8 && !collision.isInSel8())
642+
return;
643+
if (cEvtMCVtxIn10 && !collision.isVtxIn10())
644+
return;
645+
if (cEvtMCTriggerTVX && !collision.isTriggerTVX())
646+
return;
647+
if (cEvtMCRecINELgt0 && !collision.isRecINELgt0())
648+
return;
560649

561650
auto mult = collision.cent();
562651
histos.fill(HIST("Event/h1d_ft0_mult_percentile"), mult, 100);
652+
histos.fill(HIST("Event/h_ft0_vz"), collision.posZ());
563653
fillDataHistos<false, true>(tracks, tracks, mult);
564654

565655
// get MC pT-spectra
@@ -587,9 +677,9 @@ struct lambdaAnalysis_pb {
587677
}
588678

589679
for (auto const& part : resoParents) {
590-
if (abs(part.pdgCode()) != lambda1520id) // // L* pdg_code = 3124
680+
if (std::abs(part.pdgCode()) != lambda1520id) // // L* pdg_code = 3124
591681
continue;
592-
if (abs(part.y()) > 0.5) { // rapidity cut
682+
if (std::abs(part.y()) > 0.5) { // rapidity cut
593683
continue;
594684
}
595685

@@ -605,10 +695,10 @@ struct lambdaAnalysis_pb {
605695

606696
if (!pass1 || !pass2) // If we have both decay products
607697
continue;
698+
std::array<float, 3> pvec = {part.px(), part.py(), part.pz()};
699+
700+
float mass = RecoDecay::m(pvec, part.e());
608701

609-
TLorentzVector p4;
610-
p4.SetPxPyPzE(part.px(), part.py(), part.pz(), part.e());
611-
auto mass = p4.M();
612702
if (part.pdgCode() > 0)
613703
histos.fill(HIST("Analysis/h3d_gen_lstar_PM"), mass, part.pt(), mult);
614704
else

0 commit comments

Comments
 (0)