Skip to content

Commit d359453

Browse files
njacazionjacazio
authored
Refactor includes and enhance pair DCA calculations
Reorganized include statements and added new histogram registrations for filtered pairs. Updated the buildPairDCA function to streamline DCA calculations and removed redundant code.
1 parent c4f2d4e commit d359453

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)