Skip to content

Commit ebab9c6

Browse files
njacazionjacazio
authored
[ALICE3] Refactor includes and enhance pair DCA calculations (#13918)
1 parent cafae83 commit ebab9c6

File tree

1 file changed

+94
-33
lines changed

1 file changed

+94
-33
lines changed

‎ALICE3/Tasks/alice3-dilepton.cxx‎

Lines changed: 94 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -14,18 +14,21 @@
1414
/// \author s.scheid@cern.ch, daiki.sekihata@cern.ch
1515
///
1616

17-
#include "Math/Vector4D.h"
18-
// O2 includes
19-
#include "Framework/AnalysisTask.h"
20-
#include "Framework/ASoAHelpers.h"
21-
#include "Framework/runDataProcessing.h"
22-
#include "Framework/HistogramRegistry.h"
23-
#include "CommonConstants/PhysicsConstants.h"
24-
#include "Common/DataModel/TrackSelectionTables.h"
25-
#include "Framework/AnalysisDataModel.h"
26-
#include "ALICE3/DataModel/OTFTOF.h"
2717
#include "ALICE3/DataModel/OTFRICH.h"
18+
#include "ALICE3/DataModel/OTFTOF.h"
2819
#include "ALICE3/DataModel/tracksAlice3.h"
20+
#include "Common/DataModel/TrackSelectionTables.h"
21+
22+
#include <CommonConstants/PhysicsConstants.h>
23+
#include <Framework/ASoAHelpers.h>
24+
#include <Framework/AnalysisDataModel.h>
25+
#include <Framework/AnalysisTask.h>
26+
#include <Framework/HistogramRegistry.h>
27+
#include <Framework/runDataProcessing.h>
28+
29+
#include <Math/Vector4D.h>
30+
31+
#include <vector>
2932

3033
using namespace o2;
3134
using namespace o2::aod;
@@ -128,6 +131,15 @@ struct Alice3Dilepton {
128131
registry.addClone("Reconstructed/Pair/ULS/", "Reconstructed/Pair/LSpp/");
129132
registry.addClone("Reconstructed/Pair/ULS/", "Reconstructed/Pair/LSnn/");
130133

134+
registry.add("ReconstructedFiltered/Pair/ULS/Mass", "Pair Mass", kTH1F, {axisM});
135+
registry.add("ReconstructedFiltered/Pair/ULS/Pt", "Pair Pt", kTH1F, {axisPt});
136+
registry.add("ReconstructedFiltered/Pair/ULS/Eta", "Pair Eta", kTH1F, {axisEta});
137+
registry.add("ReconstructedFiltered/Pair/ULS/Phi", "Pair Phi", kTH1F, {axisPhi});
138+
registry.add("ReconstructedFiltered/Pair/ULS/Mass_Pt", "Pair Mass vs. Pt", kTH2F, {axisM, axisPt}, true);
139+
140+
registry.addClone("ReconstructedFiltered/Pair/ULS/", "ReconstructedFiltered/Pair/LSpp/");
141+
registry.addClone("ReconstructedFiltered/Pair/ULS/", "ReconstructedFiltered/Pair/LSnn/");
142+
131143
HistogramConfigSpec hs_rec{HistType::kTHnSparseF, {axisM, axisPt, axisDCAxy}, 3};
132144
registry.add("Reconstructed/Pair/ULS/hs_rec", "", hs_rec);
133145
registry.add("Reconstructed/Pair/LSpp/hs_rec", "", hs_rec);
@@ -287,6 +299,69 @@ struct Alice3Dilepton {
287299
return HFllType::kUndef;
288300
}
289301

302+
template <typename T1, typename T2>
303+
ROOT::Math::PtEtaPhiMVector buildPairDCA(T1 const& t1, T2 const& t2, float& pair_dca_xy)
304+
{
305+
306+
const float dcaXY_t1 = t1.dcaXY();
307+
const float dcaXY_t2 = t2.dcaXY();
308+
const float dcaXY_res_t1 = sqrt(t1.cYY());
309+
const float dcaXY_res_t2 = sqrt(t2.cYY());
310+
311+
pair_dca_xy = sqrt((dcaXY_t2 * dcaXY_t2 / dcaXY_res_t2 + dcaXY_t1 * dcaXY_t1 / dcaXY_res_t1) / 2.);
312+
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
313+
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
314+
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
315+
return v12;
316+
}
317+
318+
template <PairType pairtype, typename TTracks, typename TMCTracks>
319+
void FillPairRecWithPrefilter(TTracks const& tracks1, TTracks const& tracks2, TMCTracks const& /*mcParticles*/)
320+
{
321+
std::vector<uint64_t> prefilteredTracks;
322+
if constexpr (pairtype == PairType::kULS) {
323+
for (auto& [t1, t2] : combinations(soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
324+
if (!t1.has_mcParticle() || !t2.has_mcParticle()) {
325+
continue;
326+
}
327+
328+
float pair_dca_xy = 999.f;
329+
ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
330+
// prefilter for low-mass pairs
331+
if (v12.M() > 0.10) {
332+
continue;
333+
}
334+
// prefilter small opening angle pairs
335+
if (std::cos(t1.phi() - t2.phi()) < 0.99) {
336+
continue;
337+
}
338+
prefilteredTracks.push_back(t1.globalIndex());
339+
prefilteredTracks.push_back(t2.globalIndex());
340+
341+
} // end of unlike-sign pair loop
342+
343+
for (auto& [t1, t2] : combinations(soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
344+
// Skipping tracks that are in the prefiltered list
345+
if (std::find(prefilteredTracks.begin(), prefilteredTracks.end(), t1.globalIndex()) != prefilteredTracks.end()) {
346+
continue;
347+
}
348+
if (std::find(prefilteredTracks.begin(), prefilteredTracks.end(), t2.globalIndex()) != prefilteredTracks.end()) {
349+
continue;
350+
}
351+
352+
float pair_dca_xy = 999.f;
353+
ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
354+
355+
registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass"), v12.M());
356+
registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Pt"), v12.Pt());
357+
registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Eta"), v12.Eta());
358+
registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Phi"), v12.Phi() < 0.f ? v12.Phi() + TMath::TwoPi() : v12.Phi());
359+
registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt());
360+
registry.fill(HIST("ReconstructedFiltered/Pair/ULS/hs_rec"), v12.M(), v12.Pt(), pair_dca_xy);
361+
}
362+
}
363+
}
364+
290365
template <PairType pairtype, typename TTracks, typename TMCTracks>
291366
void FillPairRec(TTracks const& tracks1, TTracks const& tracks2, TMCTracks const& mcParticles)
292367
{
@@ -314,15 +389,8 @@ struct Alice3Dilepton {
314389
}
315390
// auto mother = mcparticles.iteratorAt(motherid);
316391

317-
// float dcaXY_t1 = t1.dcaXY();
318-
// float dcaXY_t2 = t2.dcaXY();
319-
// float dcaXY_res_t1 = sqrt(t1.cYY());
320-
// float dcaXY_res_t2 = sqrt(t2.cYY());
321-
322-
float pair_dca_xy = sqrt((pow(t2.dcaXY() / sqrt(t2.cYY()), 2) + pow(t1.dcaXY() / sqrt(t1.cYY()), 2)) / 2.);
323-
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
324-
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
325-
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
392+
float pair_dca_xy = 999.f;
393+
ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
326394

327395
registry.fill(HIST("Reconstructed/Pair/ULS/Mass"), v12.M());
328396
registry.fill(HIST("Reconstructed/Pair/ULS/Pt"), v12.Pt());
@@ -356,15 +424,8 @@ struct Alice3Dilepton {
356424
}
357425
// auto mother = mcparticles.iteratorAt(motherid);
358426

359-
// float dcaXY_t1 = t1.dcaXY();
360-
// float dcaXY_t2 = t2.dcaXY();
361-
// float dcaXY_res_t1 = sqrt(t1.cYY());
362-
// float dcaXY_res_t2 = sqrt(t2.cYY());
363-
364-
float pair_dca_xy = sqrt((pow(t2.dcaXY() / sqrt(t2.cYY()), 2) + pow(t1.dcaXY() / sqrt(t1.cYY()), 2)) / 2.);
365-
ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
366-
ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
367-
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
427+
float pair_dca_xy = 999.f;
428+
ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
368429

369430
if constexpr (pairtype == PairType::kLSpp) {
370431
registry.fill(HIST("Reconstructed/Pair/LSpp/Mass"), v12.M());
@@ -464,6 +525,7 @@ struct Alice3Dilepton {
464525
} // end of like-sign pair loop
465526
}
466527
}
528+
467529
// Functions for pid
468530
template <typename TTrack>
469531
bool electronIDTOF(TTrack const& track)
@@ -545,11 +607,10 @@ struct Alice3Dilepton {
545607
Partition<MyFilteredTracksMC> posTracks = o2::aod::track::signed1Pt > 0.f;
546608
Partition<MyFilteredTracksMC> negTracks = o2::aod::track::signed1Pt < 0.f;
547609

548-
void processRec(
549-
const o2::aod::Collisions& collisions,
550-
MyFilteredTracksMC const& tracks,
551-
const o2::aod::McCollisions&,
552-
const aod::McParticles& mcParticles)
610+
void processRec(const o2::aod::Collisions& collisions,
611+
MyFilteredTracksMC const& tracks,
612+
const o2::aod::McCollisions&,
613+
const aod::McParticles& mcParticles)
553614
{
554615
for (const auto& collision : collisions) {
555616
registry.fill(HIST("Reconstructed/Event/VtxX"), collision.posX());

0 commit comments

Comments
 (0)