Skip to content

[PWGDQ] added DCA analysis with global forward tracks #14981

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
iarsene merged 1 commit into from
Feb 16, 2026
+228 −57
Merged

[PWGDQ] added DCA analysis with global forward tracks #14981

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
285 changes: 228 additions & 57 deletions PWGDQ/Tasks/muonGlobalAlignment.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,8 @@ struct muonGlobalAlignment {

Configurable<bool> fEnableVertexShiftAnalysis{"cfgEnableVertexShiftAnalysis", true, "Enable the analysis of vertex shift"};
Configurable<bool> fEnableMftDcaAnalysis{"cfgEnableMftDcaAnalysis", true, "Enable the analysis of DCA-based MFT alignment"};
Configurable<bool> fEnableMftDcaExtraPlots{"cfgEnableMftDcaExtraPlots", false, "Enable additional plots for the analysis of DCA-based MFT alignment"};
Configurable<bool> fEnableGlobalFwdDcaAnalysis{"cfgEnableGlobalFwdDcaAnalysis", true, "Enable the analysis of DCA-based MFT alignment using global forward tracks"};
Configurable<bool> fEnableMftMchResidualsAnalysis{"cfgEnableMftMchResidualsAnalysis", true, "Enable the analysis of residuals between MFT tracks and MCH clusters"};

int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field
Expand Down Expand Up @@ -171,6 +173,8 @@ struct muonGlobalAlignment {
double mBzAtMftCenter{0};

HistogramRegistry registry{"registry", {}};
std::array<o2::framework::HistPtr, 10> mMftTrackEffNum;
std::array<o2::framework::HistPtr, 10> mMftTrackEffDen;

// vector of all MFT-MCH(-MID) matching candidates associated to the same MCH(-MID) track,
// to be sorted in descending order with respect to the matching quality
Expand Down Expand Up @@ -362,17 +366,22 @@ struct muonGlobalAlignment {
AxisSpec dcaxMCHAxis = {400, -10.0, 10.0, "DCA_{x} (cm)"};
AxisSpec dcayMCHAxis = {400, -10.0, 10.0, "DCA_{y} (cm)"};
AxisSpec dcazAxis = {20, -10.0, 10.0, "v_{z} (cm)"};
AxisSpec txAxis = {30, -mftLadderWidth * 15.f / 2.f, mftLadderWidth * 15.f / 2.f, "track_{x} (cm)"};
AxisSpec tyAxis = {20, -10.f, 10.f, "track_{y} (cm)"};
AxisSpec txAxis = {30 * 4, -mftLadderWidth * 15.f / 2.f, mftLadderWidth * 15.f / 2.f, "track_{x} (cm)"};
AxisSpec tyAxis = {24 * 4, -12.f, 12.f, "track_{y} (cm)"};
AxisSpec txFineAxis = {1500, -15.f, 15.f, "track_{x} (cm)"};
AxisSpec tyFineAxis = {1500, -15.f, 15.f, "track_{y} (cm)"};
AxisSpec vxAxis = {400, -0.5, 0.5, "vtx_{x} (cm)"};
AxisSpec vyAxis = {400, -0.5, 0.5, "vtx_{y} (cm)"};
AxisSpec vzAxis = {1000, -10.0, 10.0, "vtx_{z} (cm)"};
AxisSpec phiAxis = {36, -180.0, 180.0, "#phi (degrees)"};
AxisSpec momAxis = {500, 0, 100.0, "p (GeV/c)"};
AxisSpec nMftClustersAxis = {6, 5, 11, "# of clusters"};
AxisSpec mftTrackTypeAxis = {2, 0, 2, "track type"};
AxisSpec mftLayerAxis = {10, 0, 10, "layer"};
AxisSpec trackChargeSignAxis = {2, 0, 0, "sign"};
AxisSpec layersPatternAxis = {1024, 0, 1024, "layers pattern"};
AxisSpec zshiftAxis = {21, -5.25, 5.25, "z shift (mm)"};
AxisSpec chi2Axis = {500, 0, 500, "chi2"};

registry.add("vertex_y_vs_x", std::format("Vertex y vs. x").c_str(), {HistType::kTH2F, {vxAxis, vyAxis}});
registry.add("vertex_z", std::format("Vertex z").c_str(), {HistType::kTH1F, {vzAxis}});
Expand All @@ -388,12 +397,36 @@ struct muonGlobalAlignment {
}

if (fEnableMftDcaAnalysis) {
registry.add("DCA/MFT/DCA_x", "DCA(x) vs. vz, tx, ty, nclus, trackType",
HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis, mftTrackTypeAxis});
registry.add("DCA/MFT/DCA_y", "DCA(y) vs. vz, tx, ty, nclus, trackType",
HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis, mftTrackTypeAxis});
registry.add("DCA/MFT/layers", "Layers pattern vs. tx, ty, nclus, trackType",
HistType::kTHnSparseF, {layersPatternAxis, txAxis, tyAxis, nMftClustersAxis, mftTrackTypeAxis});
registry.add("DCA/MFT/DCA_x", "DCA(x) vs. vz, tx, ty, nclus",
HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis});
registry.add("DCA/MFT/DCA_y", "DCA(y) vs. vz, tx, ty, nclus",
HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis});

if (fEnableMftDcaExtraPlots) {
registry.add("DCA/MFT/layers", "Layers vs. tx, ty, nclus",
HistType::kTHnSparseF, {mftLayerAxis, txAxis, tyAxis, nMftClustersAxis});
registry.add("DCA/MFT/trackChi2", "Track #chi^{2} vs. tx, ty, nclus, layer",
HistType::kTHnSparseF, {chi2Axis, txAxis, tyAxis, nMftClustersAxis});
registry.add("DCA/MFT/trackMomentum", "Track momentum vs. tx, ty, nclus, layer",
HistType::kTHnSparseF, {momAxis, txAxis, tyAxis, nMftClustersAxis});

const int nMftLayers = 10;
for (int i = 0; i < nMftLayers; i++) {
mMftTrackEffNum[i] = registry.add((std::string("DCA/MFT/mftTrackEffNum_") + std::to_string(i)).c_str(),
(std::string("Track efficiency num. - layer ") + std::to_string(i)).c_str(),
HistType::kTH2F, {txFineAxis, tyFineAxis});
mMftTrackEffDen[i] = registry.add((std::string("DCA/MFT/mftTrackEffDen_") + std::to_string(i)).c_str(),
(std::string("Track efficiency den. - layer ") + std::to_string(i)).c_str(),
HistType::kTH2F, {txFineAxis, tyFineAxis});
}
}
}

if (fEnableGlobalFwdDcaAnalysis) {
registry.add("DCA/GlobalFwd/DCA_x", "DCA(x) vs. vz, tx, ty, nclus",
HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis});
registry.add("DCA/GlobalFwd/DCA_y", "DCA(y) vs. vz, tx, ty, nclus",
HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis});
}

if (fEnableMftMchResidualsAnalysis) {
Expand Down Expand Up @@ -934,6 +967,68 @@ struct muonGlobalAlignment {
fwdtrack.setCovariances(transformedTrack.getCovariances());
}

template <typename T>
T UpdateTrackMomentum(const T& track, const double p, int sign)
{
double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());
double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));

SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt};
std::vector<double> v1{0, 0, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, 0, 0, 0};
SMatrix55 tcovs(v1.begin(), v1.end());

T newTrack;
newTrack.setParameters(tpars);
newTrack.setZ(track.z());
newTrack.setCovariances(tcovs);

return newTrack;
}

template <typename T>
T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom)
{
double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());
double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
double sign = track4mom.getCharge();

SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt};
std::vector<double> v1{0, 0, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, 0, 0, 0};
SMatrix55 tcovs(v1.begin(), v1.end());

T newTrack;
newTrack.setParameters(tpars);
newTrack.setZ(track.z());
newTrack.setCovariances(tcovs);

return track;
}

void UpdateTrackMomentum(o2::mch::TrackParam& track, const o2::mch::TrackParam& track4mom)
{
double pRatio = track.p() / track4mom.p();
double newInvBendMom = track.getInverseBendingMomentum() * pRatio;
track.setInverseBendingMomentum(newInvBendMom);
track.setCharge(track4mom.getCharge());
}

void UpdateTrackMomentum(o2::track::TrackParCovFwd& track, const o2::mch::TrackParam& track4mom)
{
auto newTrackMCH = FwdtoMCH(track);
UpdateTrackMomentum(newTrackMCH, track4mom);
auto newTrack = MCHtoFwd(newTrackMCH);

track.setParameters(newTrack.getParameters());
track.setZ(newTrack.getZ());
track.setCovariances(newTrack.getCovariances());
}

o2::dataformats::GlobalFwdTrack PropagateMCHParam(mch::TrackParam mchTrack, const double z)
{
float absFront = -90.f;
Expand Down Expand Up @@ -1041,7 +1136,7 @@ struct muonGlobalAlignment {
}

template <class TMFT, class C>
o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const C& collision, float zshift = 0)
o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const C& collision, float zshift)
{
static double Bz = -10001;
double chi2 = mftTrack.chi2();
Expand Down Expand Up @@ -1084,55 +1179,53 @@ struct muonGlobalAlignment {
return propmuon;
}

template <typename T>
T UpdateTrackMomentum(const T& track, const double p, int sign)
template <class TMFT, class TMUON, class C>
o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const TMUON& mchTrack, const C& collision, float zshift)
{
double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());
double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));

SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt};
static double Bz = -10001;
double chi2 = mftTrack.chi2();
double phiCorrDeg = 0;
double phiCorr = phiCorrDeg * TMath::Pi() / 180.f;
double tR = std::hypot(mftTrack.x(), mftTrack.y());
double tphi = std::atan2(mftTrack.y(), mftTrack.x());
double tx = std::cos(tphi + phiCorr) * tR;
double ty = std::sin(tphi + phiCorr) * tR;
SMatrix5 tpars = {tx, ty, mftTrack.phi() + phiCorr, mftTrack.tgl(), mftTrack.signed1Pt()};
std::vector<double> v1{0, 0, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, 0, 0, 0};
SMatrix55 tcovs(v1.begin(), v1.end());
o2::track::TrackParCovFwd fwdtrack{mftTrack.z(), tpars, tcovs, chi2};
if (configMFTAlignmentCorrections.fEnableMFTAlignmentCorrections) {
TransformMFT(fwdtrack);
}

T newTrack;
newTrack.setParameters(tpars);
newTrack.setZ(track.z());
newTrack.setCovariances(tcovs);

return newTrack;
}

template <typename T>
T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom)
{
double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());
double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
double sign = track4mom.getCharge();
// extrapolation with MCH tools
auto mchTrackAtMFT = FwdtoMCH(FwdToTrackPar(mchTrack));
o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z());
UpdateTrackMomentum(fwdtrack, mchTrackAtMFT);

SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt};
std::vector<double> v1{0, 0, 0, 0, 0,
0, 0, 0, 0, 0,
0, 0, 0, 0, 0};
SMatrix55 tcovs(v1.begin(), v1.end());
// double propVec[3] = {};
// propVec[0] = collision.posX() - mftTrack.x();
// propVec[1] = collision.posY() - mftTrack.y();
// propVec[2] = collision.posZ() - mftTrack.z();

T newTrack;
newTrack.setParameters(tpars);
newTrack.setZ(track.z());
newTrack.setCovariances(tcovs);
// double centerZ[3] = {mftTrack.x() + propVec[0] / 2.,
// mftTrack.y() + propVec[1] / 2.,
// mftTrack.z() + propVec[2] / 2.};
if (Bz < -10000) {
double centerZ[3] = {0, 0, -45.f / 2.f};
o2::field::MagneticField* field = static_cast<o2::field::MagneticField*>(TGeoGlobalMagField::Instance()->GetField());
Bz = field->getBz(centerZ);
}
fwdtrack.propagateToZ(collision.posZ() - zshift, Bz);

return track;
}
o2::dataformats::GlobalFwdTrack propmuon;
propmuon.setParameters(fwdtrack.getParameters());
propmuon.setZ(fwdtrack.getZ());
propmuon.setCovariances(fwdtrack.getCovariances());

void UpdateTrackMomentum(o2::mch::TrackParam& track, const o2::mch::TrackParam& track4mom)
{
double pRatio = track.p() / track4mom.p();
double newInvBendMom = track.getInverseBendingMomentum() * pRatio;
track.setInverseBendingMomentum(newInvBendMom);
track.setCharge(track4mom.getCharge());
return propmuon;
}

template <typename TMCH, typename TMFT>
Expand Down Expand Up @@ -1170,7 +1263,7 @@ struct muonGlobalAlignment {

void FillDCAPlots(MyEvents const& collisions,
MyBCs const& bcs,
MyMuonsWithCov const& /*muonTracks*/,
MyMuonsWithCov const& muonTracks,
MyMFTs const& mftTracks,
const std::map<uint64_t, CollisionInfo>& collisionInfos)
{
Expand Down Expand Up @@ -1203,6 +1296,10 @@ struct muonGlobalAlignment {
for (auto mftIndex : mftTrackIds) {
auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex);

if (mftTrack.isCA()) {
continue;
}

bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5);
if (!isGoodMFT)
continue;
Expand All @@ -1212,26 +1309,55 @@ struct muonGlobalAlignment {
double dcay = mftTrackAtDCA.getY() - collision.posY();
double phi = mftTrack.phi() * 180 / TMath::Pi();
int mftNclusters = mftTrack.nClusters();
int mftTrackType = mftTrack.isCA() ? 1 : 0;
double chi2NDF = static_cast<double>(mftNclusters) * 2 - 5;

const int nMftLayers = 10;
int layerPattern = 0;
std::array<bool, 10> firedLayers;
for (int layer = 0; layer < nMftLayers; layer++) {
if ((mftTrack.mftClusterSizesAndTrackFlags() >> (layer * 6)) & 0x3F) {
layerPattern += (1 << layer);
firedLayers[layer] = true;
} else {
firedLayers[layer] = false;
}
}

if (fEnableMftDcaAnalysis) {
registry.get<TH2>(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay);
registry.get<THnSparse>(HIST("DCA/MFT/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters, mftTrackType);
registry.get<THnSparse>(HIST("DCA/MFT/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters, mftTrackType);
const int nMftLayers = 10;
if (fEnableMftDcaExtraPlots) {
for (int i = 0; i < nMftLayers; i++) {
if (firedLayers[i]) {
registry.get<THnSparse>(HIST("DCA/MFT/trackChi2"))->Fill(mftTrack.chi2() / chi2NDF, mftTrack.x(), mftTrack.y(), mftNclusters, i);
}
}
}

registry.get<THnSparse>(HIST("DCA/MFT/layers"))->Fill(layerPattern, mftTrack.x(), mftTrack.y(), mftNclusters, mftTrackType);
if (mftTrack.chi2() <= fTrackChi2MftUp) {
registry.get<TH2>(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay);
registry.get<THnSparse>(HIST("DCA/MFT/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters);
registry.get<THnSparse>(HIST("DCA/MFT/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters);

if (fEnableMftDcaExtraPlots) {
if (mftNclusters >= 6) {
for (int i = 0; i < nMftLayers; i++) {
auto mftTrackAtLayer = PropagateMFT(mftTrack, o2::mft::constants::mft::LayerZCoordinate()[i]);
std::get<std::shared_ptr<TH2>>(mMftTrackEffDen[i])->Fill(mftTrackAtLayer.getX(), mftTrackAtLayer.getY());
if (firedLayers[i]) {
std::get<std::shared_ptr<TH2>>(mMftTrackEffNum[i])->Fill(mftTrackAtLayer.getX(), mftTrackAtLayer.getY());
registry.get<THnSparse>(HIST("DCA/MFT/layers"))->Fill(i, mftTrack.x(), mftTrack.y(), mftNclusters);
}
}
}
for (int i = 0; i < nMftLayers; i++) {
if (firedLayers[i]) {
registry.get<THnSparse>(HIST("DCA/MFT/trackMomentum"))->Fill(mftTrack.p(), mftTrack.x(), mftTrack.y(), mftNclusters, i);
}
}
}
}
}

if (fEnableVertexShiftAnalysis) {
if (std::fabs(collision.posZ()) < 1.f) {
if (mftTrack.chi2() <= fTrackChi2MftUp && std::fabs(collision.posZ()) < 1.f) {
float zshift[21] = {// in millimeters
-5.0, -4.5, -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0,
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0};
Expand All @@ -1246,6 +1372,51 @@ struct muonGlobalAlignment {
}
}
}

if (fEnableGlobalFwdDcaAnalysis) {
// loop over global muon tracks
for (auto& [muonIndex, globalTracksVector] : collisionInfo.globalMuonTracks) {
auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0]);
const auto& mchTrack = muonTrack.template matchMCHTrack_as<MyMuonsWithCov>();
const auto& mftTrack = muonTrack.template matchMFTTrack_as<MyMFTs>();

if (muonTrack.chi2MatchMCHMFT() < 50) {
continue;
}

if (globalTracksVector.size() > 1) {
auto const& muonTrack2 = muonTracks.rawIteratorAt(globalTracksVector[1]);
double dchi2 = muonTrack2.chi2MatchMCHMFT() - muonTrack.chi2MatchMCHMFT();

if (dchi2 < 50) {
continue;
}
}

if (mftTrack.isCA()) {
continue;
}

bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5);
if (!isGoodMFT) {
continue;
}

bool isGoodMuon = IsGoodMuon(mchTrack, collision, fTrackChi2MchUp, 0.f, fPtMchLow, {fEtaMftLow, fEtaMftUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp);
if (!isGoodMuon)
continue;

auto mftTrackAtDCA = PropagateMFTToDCA(mftTrack, mchTrack, collision, fVertexZshift);
double dcax = mftTrackAtDCA.getX() - collision.posX();
double dcay = mftTrackAtDCA.getY() - collision.posY();
int mftNclusters = mftTrack.nClusters();

if (mftTrack.chi2() <= fTrackChi2MftUp) {
registry.get<THnSparse>(HIST("DCA/GlobalFwd/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters);
registry.get<THnSparse>(HIST("DCA/GlobalFwd/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters);
}
}
}
}
}

Expand Down
Loading