Skip to content

Commit 66ee41c

Browse files
rolavickrolavick
committed
new propagation to vertex for global tracks
1 parent 09779f0 commit 66ee41c

File tree

1 file changed

+119
-12
lines changed

1 file changed

+119
-12
lines changed

‎PWGUD/TableProducer/upcCandProducerGlobalMuon.cxx‎

Lines changed: 119 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,8 @@
2525
#include "Framework/AnalysisDataModel.h"
2626
#include "Framework/AnalysisTask.h"
2727
#include "Framework/runDataProcessing.h"
28+
#include "Common/Core/fwdtrackUtilities.h"
29+
#include "CommonConstants/GeomConstants.h"
2830
#include "GlobalTracking/MatchGlobalFwd.h"
2931
#include "MCHTracking/TrackExtrap.h"
3032
#include "ReconstructionDataFormats/TrackFwd.h"
@@ -71,6 +73,13 @@ struct UpcCandProducerGlobalMuon {
7173
Configurable<float> fMaxEtaMFT{"fMaxEtaMFT", -2.5, "Maximum eta for MFT acceptance"};
7274
Configurable<bool> fSaveMFTClusters{"fSaveMFTClusters", true, "Save MFT cluster information"};
7375

76+
// Ambiguous track propagation configurables
77+
Configurable<bool> fApplyZShiftFromCCDB{"fApplyZShiftFromCCDB", false, "Apply z-shift from CCDB for global muon propagation"};
78+
Configurable<std::string> fZShiftPath{"fZShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z-shift"};
79+
Configurable<float> fManualZShift{"fManualZShift", 0.0f, "Manual z-shift for global muon propagation to PV"};
80+
Configurable<float> fMaxDCAxy{"fMaxDCAxy", 999.f, "Maximum DCAxy for global muon track selection (cm)"};
81+
Configurable<int> fBcWindowCollision{"fBcWindowCollision", 4, "BC window for collision search for DCA-based vertex assignment"};
82+
7483
using ForwardTracks = o2::soa::Join<o2::aod::FwdTracks, o2::aod::FwdTracksCov>;
7584

7685
HistogramRegistry histRegistry{"HistRegistry", {}, OutputObjHandlingPolicy::AnalysisObject};
@@ -80,10 +89,15 @@ struct UpcCandProducerGlobalMuon {
8089
o2::ccdb::CcdbApi fCCDBApi;
8190
o2::globaltracking::MatchGlobalFwd fMatching;
8291

92+
// Ambiguous track propagation members
93+
float fBz{0}; // Magnetic field at MFT center
94+
static constexpr double fCcenterMFT[3] = {0, 0, -61.4}; // Field evaluation point at center of MFT
95+
float fZShift{0}; // z-vertex shift for forward track propagation
96+
8397
void init(InitContext&)
8498
{
8599
fUpcCuts = (UPCCutparHolder)fUpcCutsConf;
86-
if (fUpcCuts.getUseFwdCuts()) {
100+
if (fUpcCuts.getUseFwdCuts() || fEnableMFT) {
87101
fCCDB->setURL("http://alice-ccdb.cern.ch");
88102
fCCDB->setCaching(true);
89103
fCCDB->setLocalObjectValidityChecking();
@@ -110,6 +124,12 @@ struct UpcCandProducerGlobalMuon {
110124
const AxisSpec axisEta{100, -4.0, -2.0, "#eta"};
111125
histRegistry.add("hEtaMFT", "MFT track eta", kTH1F, {axisEta});
112126
histRegistry.add("hEtaGlobal", "Global track eta", kTH1F, {axisEta});
127+
128+
const AxisSpec axisDCAxy{200, 0., 10., "DCA_{xy} (cm)"};
129+
const AxisSpec axisDCAz{200, -10., 10., "DCA_{z} (cm)"};
130+
histRegistry.add("hDCAxyGlobal", "DCAxy of global tracks to best collision", kTH1F, {axisDCAxy});
131+
histRegistry.add("hDCAzGlobal", "DCAz of global tracks to best collision", kTH1F, {axisDCAz});
132+
histRegistry.add("hNCompatColls", "Number of compatible collisions per global track", kTH1F, {{21, -0.5, 20.5}});
113133
}
114134
}
115135

@@ -436,6 +456,19 @@ struct UpcCandProducerGlobalMuon {
436456
return (eta > fMinEtaMFT && eta < fMaxEtaMFT);
437457
}
438458

459+
// Propagate global muon track to collision vertex using helix propagation
460+
// and compute DCA (adapted from ambiguousTrackPropagation)
461+
// Returns {DCAxy, DCAz, DCAx, DCAy}
462+
std::array<double, 4> propagateGlobalToDCA(ForwardTracks::iterator const& track,
463+
double colX, double colY, double colZ)
464+
{
465+
o2::track::TrackParCovFwd trackPar = o2::aod::fwdtrackutils::getTrackParCovFwdShift(track, fZShift);
466+
std::array<double, 3> dcaOrig;
467+
trackPar.propagateToDCAhelix(fBz, {colX, colY, colZ}, dcaOrig);
468+
double dcaXY = std::sqrt(dcaOrig[0] * dcaOrig[0] + dcaOrig[1] * dcaOrig[1]);
469+
return {dcaXY, dcaOrig[2], dcaOrig[0], dcaOrig[1]};
470+
}
471+
439472
void createCandidates(ForwardTracks const& fwdTracks,
440473
o2::aod::FwdTrkCls const& fwdTrkCls,
441474
o2::aod::AmbiguousFwdTracks const& ambFwdTracks,
@@ -451,7 +484,7 @@ struct UpcCandProducerGlobalMuon {
451484
using o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack;
452485

453486
int32_t runNumber = bcs.iteratorAt(0).runNumber();
454-
if (fUpcCuts.getUseFwdCuts()) {
487+
if (fUpcCuts.getUseFwdCuts() || fEnableMFT) {
455488
if (runNumber != fRun) {
456489
fRun = runNumber;
457490
std::map<std::string, std::string> metadata;
@@ -461,7 +494,28 @@ struct UpcCandProducerGlobalMuon {
461494
o2::base::Propagator::initFieldFromGRP(grpmag);
462495
if (!o2::base::GeometryManager::isGeometryLoaded())
463496
fCCDB->get<TGeoManager>("GLO/Config/GeometryAligned");
464-
o2::mch::TrackExtrap::setField();
497+
if (fUpcCuts.getUseFwdCuts()) {
498+
o2::mch::TrackExtrap::setField();
499+
}
500+
// Initialize MFT magnetic field and z-shift for ambiguous track propagation
501+
if (fEnableMFT) {
502+
o2::field::MagneticField* field = static_cast<o2::field::MagneticField*>(TGeoGlobalMagField::Instance()->GetField());
503+
fBz = field->getBz(fCcenterMFT);
504+
LOG(info) << "Magnetic field at MFT center: bZ = " << fBz;
505+
if (fApplyZShiftFromCCDB) {
506+
auto* zShift = fCCDB->getForTimeStamp<std::vector<float>>(fZShiftPath, ts);
507+
if (zShift != nullptr && !zShift->empty()) {
508+
fZShift = (*zShift)[0];
509+
LOGF(info, "z-shift from CCDB: %f cm", fZShift);
510+
} else {
511+
fZShift = 0;
512+
LOGF(info, "z-shift not found in CCDB path %s, set to 0 cm", fZShiftPath.value);
513+
}
514+
} else {
515+
fZShift = fManualZShift;
516+
LOGF(info, "z-shift manually set to %f cm", fZShift);
517+
}
518+
}
465519
}
466520
}
467521

@@ -533,6 +587,15 @@ struct UpcCandProducerGlobalMuon {
533587
}
534588
}
535589

590+
// Map global BC to collisions for DCA-based vertex assignment
591+
std::map<uint64_t, std::vector<int64_t>> mapGlobalBCtoCollisions;
592+
if (fEnableMFT) {
593+
for (const auto& col : collisions) {
594+
uint64_t gbc = vGlobalBCs[col.bcId()];
595+
mapGlobalBCtoCollisions[gbc].push_back(col.globalIndex());
596+
}
597+
}
598+
536599
std::vector<int64_t> selTrackIds{}; // NEW: For cluster saving
537600

538601
int32_t candId = 0;
@@ -552,34 +615,77 @@ struct UpcCandProducerGlobalMuon {
552615
continue;
553616
}
554617

555-
uint16_t numContrib = 0;
556618
auto& vGlobalIds = gbc_globalids.second;
557619

558-
// Check if we have global tracks (with MFT)
620+
// Step 1: Find best collision vertex using DCA-based propagation
621+
// (adapted from ambiguousTrackPropagation processMFTReassoc3D)
622+
float bestVtxX = 0., bestVtxY = 0., bestVtxZ = 0.;
623+
double bestAvgDCA = 999.;
624+
bool hasVertex = false;
625+
int nCompatColls = 0;
626+
627+
for (int dbc = -static_cast<int>(fBcWindowCollision); dbc <= static_cast<int>(fBcWindowCollision); dbc++) {
628+
uint64_t searchBC = globalBcGlobal + dbc;
629+
auto itCol = mapGlobalBCtoCollisions.find(searchBC);
630+
if (itCol == mapGlobalBCtoCollisions.end())
631+
continue;
632+
for (auto colIdx : itCol->second) {
633+
nCompatColls++;
634+
const auto& col = collisions.iteratorAt(colIdx);
635+
double sumDCAxy = 0.;
636+
int nTracks = 0;
637+
for (const auto& iglobal : vGlobalIds) {
638+
const auto& trk = fwdTracks.iteratorAt(iglobal);
639+
auto dca = propagateGlobalToDCA(trk, col.posX(), col.posY(), col.posZ());
640+
sumDCAxy += dca[0];
641+
nTracks++;
642+
}
643+
double avgDCA = nTracks > 0 ? sumDCAxy / nTracks : 999.;
644+
if (!hasVertex || avgDCA < bestAvgDCA) {
645+
bestAvgDCA = avgDCA;
646+
bestVtxX = col.posX();
647+
bestVtxY = col.posY();
648+
bestVtxZ = col.posZ();
649+
hasVertex = true;
650+
}
651+
}
652+
}
653+
654+
histRegistry.fill(HIST("hNCompatColls"), nCompatColls);
655+
656+
// Step 2: Select tracks with DCA quality cut
559657
std::vector<int64_t> tracksToSave;
560658
for (const auto& iglobal : vGlobalIds) {
561659
const auto& trk = fwdTracks.iteratorAt(iglobal);
562660

661+
// Apply DCA cut using best collision vertex
662+
if (hasVertex) {
663+
auto dca = propagateGlobalToDCA(trk, bestVtxX, bestVtxY, bestVtxZ);
664+
histRegistry.fill(HIST("hDCAxyGlobal"), dca[0]);
665+
histRegistry.fill(HIST("hDCAzGlobal"), dca[1]);
666+
if (dca[0] > static_cast<double>(fMaxDCAxy))
667+
continue;
668+
}
669+
563670
// Check MFT acceptance and decide which track to use
564671
if (isInMFTAcceptance(trk.eta())) {
565672
// Inside MFT acceptance - use global track
566673
tracksToSave.push_back(iglobal);
567674
histRegistry.fill(HIST("hEtaMFT"), trk.eta());
568675
} else {
569676
// Outside MFT acceptance - look for MCH-MID counterpart
570-
// Find the corresponding MCH-MID track at the same BC
571677
auto itMid = mapGlobalBcsWithMCHMIDTrackIds.find(globalBcGlobal);
572678
if (itMid != mapGlobalBcsWithMCHMIDTrackIds.end()) {
573-
// Use MCH-MID track instead
574679
if (!itMid->second.empty()) {
575680
tracksToSave.push_back(itMid->second[0]);
576-
itMid->second.erase(itMid->second.begin()); // Remove used track
681+
itMid->second.erase(itMid->second.begin());
577682
}
578683
}
579684
}
580685
}
581686

582-
// Write selected tracks
687+
// Step 3: Write tracks and event candidate with actual vertex position
688+
uint16_t numContrib = 0;
583689
for (const auto& trkId : tracksToSave) {
584690
if (!addToFwdTable(candId, trkId, globalBcGlobal, 0., fwdTracks, mcFwdTrackLabels))
585691
continue;
@@ -590,7 +696,7 @@ struct UpcCandProducerGlobalMuon {
590696
if (numContrib < 1)
591697
continue;
592698

593-
eventCandidates(globalBcGlobal, runNumber, 0., 0., 0., 0, numContrib, 0, 0);
699+
eventCandidates(globalBcGlobal, runNumber, bestVtxX, bestVtxY, bestVtxZ, 0, numContrib, 0, 0);
594700
std::vector<float> amplitudesV0A{};
595701
std::vector<int8_t> relBCsV0A{};
596702
std::vector<float> amplitudesT0A{};
@@ -683,8 +789,9 @@ struct UpcCandProducerGlobalMuon {
683789
vAmbFwdTrackIndexBCs.clear();
684790
mapGlobalBcsWithMCHMIDTrackIds.clear();
685791
mapGlobalBcsWithMCHTrackIds.clear();
686-
mapGlobalBcsWithGlobalTrackIds.clear(); // NEW
687-
selTrackIds.clear(); // NEW
792+
mapGlobalBcsWithGlobalTrackIds.clear();
793+
mapGlobalBCtoCollisions.clear();
794+
selTrackIds.clear();
688795
}
689796

690797
void processFwd(ForwardTracks const& fwdTracks,

0 commit comments

Comments
 (0)