Skip to content

Commit b4f8e55

Browse files
mcoquet642mcoquet642
authored
[PWGDQ] Code cleanup and logic improvements for muon processing in tableMaker(MC)_withAssoc (#17893)
1 parent e76d66c commit b4f8e55

3 files changed

Lines changed: 119 additions & 121 deletions

File tree

‎PWGDQ/Core/VarManager.h‎

Lines changed: 39 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -1449,8 +1449,6 @@ class VarManager : public TObject
14491449
template <typename T, typename C>
14501450
static o2::track::TrackParCovFwd PropagateFwd(const T& track, const C& cov, float z);
14511451
template <uint32_t fillMap, typename T, typename C>
1452-
static void FillMuonPDca(const T& muon, const C& collision, float* values = nullptr);
1453-
template <uint32_t fillMap, typename T, typename C>
14541452
static void FillPropagateMuon(const T& muon, const C& collision, float* values = nullptr);
14551453
template <typename T>
14561454
static void FillBC(T const& bc, float* values = nullptr);
@@ -1837,11 +1835,7 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C
18371835
o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwd3DShift(muon, xShift, yShift, zShift, muon);
18381836
o2::dataformats::GlobalFwdTrack propmuon;
18391837
if (static_cast<int>(muon.trackType()) > 2) {
1840-
o2::dataformats::GlobalFwdTrack track;
1841-
track.setParameters(fwdtrack.getParameters());
1842-
track.setZ(fwdtrack.getZ());
1843-
track.setCovariances(fwdtrack.getCovariances());
1844-
auto mchTrack = mMatching.FwdtoMCH(track);
1838+
auto mchTrack = mMatching.FwdtoMCH(fwdtrack);
18451839

18461840
if (endPoint == kToVertex) {
18471841
o2::mch::TrackExtrap::extrapToVertex(mchTrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY());
@@ -1856,17 +1850,12 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C
18561850
o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, fgzMatching);
18571851
}
18581852

1859-
auto proptrack = mMatching.MCHtoFwd(mchTrack);
1860-
propmuon.setParameters(proptrack.getParameters());
1861-
propmuon.setZ(proptrack.getZ());
1862-
propmuon.setCovariances(proptrack.getCovariances());
1853+
propmuon = mMatching.MCHtoFwd(mchTrack);
18631854

18641855
} else if (static_cast<int>(muon.trackType()) < 2) {
18651856
std::array<double, 3> dcaInfOrig{999.f, 999.f, 999.f};
18661857
fwdtrack.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dcaInfOrig);
1867-
propmuon.setParameters(fwdtrack.getParameters());
1868-
propmuon.setZ(fwdtrack.getZ());
1869-
propmuon.setCovariances(fwdtrack.getCovariances());
1858+
propmuon = fwdtrack;
18701859
}
18711860
return propmuon;
18721861
}
@@ -1879,25 +1868,6 @@ o2::track::TrackParCovFwd VarManager::PropagateFwd(const T& track, const C& cov,
18791868
return fwdtrack;
18801869
}
18811870

1882-
template <uint32_t fillMap, typename T, typename C>
1883-
void VarManager::FillMuonPDca(const T& muon, const C& collision, float* values)
1884-
{
1885-
if (!values) {
1886-
values = fgValues;
1887-
}
1888-
1889-
if constexpr ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0) {
1890-
1891-
o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(muon, collision);
1892-
o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(muon, collision, kToDCA);
1893-
1894-
float dcaX = (propmuonAtDCA.getX() - collision.posX());
1895-
float dcaY = (propmuonAtDCA.getY() - collision.posY());
1896-
float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY);
1897-
values[kMuonPDca] = muon.p() * dcaXY;
1898-
}
1899-
}
1900-
19011871
template <uint32_t fillMap, typename T, typename C>
19021872
void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* values)
19031873
{
@@ -1921,20 +1891,6 @@ void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* val
19211891
values[kTgl] = propmuon.getTgl();
19221892
values[kPhi] = propmuon.getPhi();
19231893

1924-
// Redo propagation only for muon tracks
1925-
// propagation of MFT tracks alredy done in fwdtrack-extention task
1926-
if (static_cast<int>(muon.trackType()) > 2) {
1927-
o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(muon, collision, kToDCA);
1928-
o2::dataformats::GlobalFwdTrack propmuonAtRabs = PropagateMuon(muon, collision, kToRabs);
1929-
float dcaX = (propmuonAtDCA.getX() - collision.posX());
1930-
float dcaY = (propmuonAtDCA.getY() - collision.posY());
1931-
values[kMuonDCAx] = dcaX;
1932-
values[kMuonDCAy] = dcaY;
1933-
double xAbs = propmuonAtRabs.getX();
1934-
double yAbs = propmuonAtRabs.getY();
1935-
values[kMuonRAtAbsorberEnd] = std::sqrt(xAbs * xAbs + yAbs * yAbs);
1936-
}
1937-
19381894
const SMatrix55& cov = propmuon.getCovariances();
19391895
values[kMuonCXX] = cov(0, 0);
19401896
values[kMuonCXY] = cov(1, 0);
@@ -1971,6 +1927,7 @@ void VarManager::FillGlobalMuonRefit(T1 const& muontrack, T2 const& mfttrack, co
19711927
double pz = propmuon.getP() * std::cos(o2::constants::math::PIHalf - std::atan(mfttrack.tgl()));
19721928
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
19731929
auto mftprop = o2::aod::fwdtrackutils::getTrackParCovFwd3DShift(mfttrack, xShift, yShift, zShift);
1930+
mftprop.setInvQPt(static_cast<double>(muontrack.sign()) / pt);
19741931
values[kX] = mftprop.getX();
19751932
values[kY] = mftprop.getY();
19761933
values[kZ] = mftprop.getZ();
@@ -1979,6 +1936,12 @@ void VarManager::FillGlobalMuonRefit(T1 const& muontrack, T2 const& mfttrack, co
19791936
values[kPz] = pz;
19801937
values[kEta] = mftprop.getEta();
19811938
values[kPhi] = mftprop.getPhi();
1939+
1940+
// Helix DCA of the refitted global track w.r.t. the associated collision
1941+
std::array<double, 3> dca{999., 999., 999.};
1942+
mftprop.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dca);
1943+
values[kMuonDCAx] = static_cast<float>(dca[0]);
1944+
values[kMuonDCAy] = static_cast<float>(dca[1]);
19821945
}
19831946
}
19841947

@@ -2006,6 +1969,12 @@ void VarManager::FillGlobalMuonRefitCov(T1 const& muontrack, T2 const& mfttrack,
20061969
values[kPz] = globalRefit.getPz();
20071970
values[kEta] = globalRefit.getEta();
20081971
values[kPhi] = globalRefit.getPhi();
1972+
1973+
// Helix DCA of the covariance-refitted global track w.r.t. the associated collision
1974+
std::array<double, 3> dca{999., 999., 999.};
1975+
globalRefit.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dca);
1976+
values[kMuonDCAx] = static_cast<float>(dca[0]);
1977+
values[kMuonDCAy] = static_cast<float>(dca[1]);
20091978
}
20101979
}
20111980
}
@@ -3477,13 +3446,30 @@ void VarManager::FillTrackCollision(T const& track, C const& collision, float* v
34773446
}
34783447
}
34793448
if constexpr ((fillMap & MuonCov) > 0 || (fillMap & MuonCovRealign) > 0 || (fillMap & ReducedMuonCov) > 0) {
3449+
float dcaX = 999.f;
3450+
float dcaY = 999.f;
3451+
if (static_cast<int>(track.trackType()) <= 2) {
3452+
// Global / MCH-MID: helix DCA only (for globals, kMuonPDca is filled from the matched MCH in skimMuons)
3453+
float xShift = 0.f;
3454+
float yShift = 0.f;
3455+
float zShift = 0.f;
3456+
GetFwdShiftForY(track.y(), xShift, yShift, zShift);
3457+
o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwd3DShift(track, xShift, yShift, zShift, track);
3458+
std::array<double, 3> dca{999., 999., 999.};
3459+
fwdtrack.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dca);
3460+
dcaX = static_cast<float>(dca[0]);
3461+
dcaY = static_cast<float>(dca[1]);
3462+
} else {
3463+
// MCH standalone: DCA and pDCA from MCH extrapolation
3464+
o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(track, collision, kToDCA);
3465+
dcaX = propmuonAtDCA.getX() - collision.posX();
3466+
dcaY = propmuonAtDCA.getY() - collision.posY();
3467+
float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY);
3468+
values[kMuonPDca] = track.p() * dcaXY;
3469+
}
34803470

3481-
o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(track, collision, kToDCA);
3482-
3483-
float dcaX = (propmuonAtDCA.getX() - collision.posX());
3484-
float dcaY = (propmuonAtDCA.getY() - collision.posY());
3485-
float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY);
3486-
values[kMuonPDca] = track.p() * dcaXY;
3471+
values[kMuonDCAx] = dcaX;
3472+
values[kMuonDCAy] = dcaY;
34873473
}
34883474
}
34893475

‎PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx‎

Lines changed: 41 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,7 @@ template <typename TMap>
144144
void PrintBitMap(TMap map, int nbits)
145145
{
146146
for (int i = 0; i < nbits; i++) {
147-
cout << ((map & (TMap(1) << i)) > 0 ? "1" : "0");
147+
LOG(info) << ((map & (TMap(1) << i)) > 0 ? "1" : "0");
148148
}
149149
}
150150
*/
@@ -234,9 +234,6 @@ struct TableMakerMC {
234234
Configurable<std::string> fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
235235
Configurable<std::string> fGeoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
236236
Configurable<std::string> fGrpMagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
237-
Configurable<std::string> fZShiftPath{"zShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"};
238-
Configurable<bool> fUseRemoteZShift{"cfgUseRemoteZShift", false, "Enable getting Zshift from ccdb"};
239-
Configurable<float> fManualZShift{"cfgManualZShift", 0.f, "Manual value for the Zshift for muons."};
240237
Configurable<std::string> fGrpMagPathRun2{"grpmagPathRun2", "GLO/GRP/GRP", "CCDB path of the GRPObject (Usage for Run 2)"};
241238
Configurable<int64_t> timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"};
242239
} fConfigCCDB;
@@ -580,17 +577,16 @@ struct TableMakerMC {
580577
/*if ((std::abs(mctrack.pdgCode())>400 && std::abs(mctrack.pdgCode())<599) ||
581578
(std::abs(mctrack.pdgCode())>4000 && std::abs(mctrack.pdgCode())<5999) ||
582579
(mcflags > 0)) {
583-
cout << ">>>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : "
584-
<< mctrack.globalIndex() << " / " << mctrack.pdgCode() << " / "
585-
<< mctrack.getProcess() << " / " << mctrack.getGenStatusCode() << " / " << mctrack.getHepMCStatusCode() << " / " << mctrack.isPhysicalPrimary() << endl;
586-
cout << ">>>>>>>>>>>>>>>>>>>>>>> track bitmap: ";
580+
LOG(info) << ">>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : "
581+
<< mctrack.globalIndex() << " / " << mctrack.pdgCode() << " / "
582+
<< mctrack.getProcess() << " / " << mctrack.getGenStatusCode() << " / " << mctrack.getHepMCStatusCode() << " / " << mctrack.isPhysicalPrimary();
583+
LOG(info) << ">>>>>>>>>>>>>>>>>>>>>> track bitmap: ";
587584
PrintBitMap(mcflags, 16);
588-
cout << endl;
589585
if (mctrack.has_mothers()) {
590586
for (const auto& m : mctrack.mothersIds()) {
591587
if (m < mcTracks.size()) { // protect against bad mother indices
592588
auto aMother = mcTracks.rawIteratorAt(m);
593-
cout << "<<<<<< mother idx / pdg: " << m << " / " << aMother.pdgCode() << endl;
589+
LOG(info) << "<<<<<< mother idx / pdg: " << m << " / " << aMother.pdgCode();
594590
}
595591
}
596592
}
@@ -600,7 +596,7 @@ struct TableMakerMC {
600596
601597
if (d < mcTracks.size()) { // protect against bad daughter indices
602598
auto aDaughter = mcTracks.rawIteratorAt(d);
603-
cout << "<<<<<< daughter idx / pdg: " << d << " / " << aDaughter.pdgCode() << endl;
599+
LOG(info) << "<<<<<< daughter idx / pdg: " << d << " / " << aDaughter.pdgCode();
604600
}
605601
}
606602
}
@@ -1133,22 +1129,29 @@ struct TableMakerMC {
11331129
VarManager::FillTrack<TMuonFillMap>(muon);
11341130
// NOTE: If a muon is associated to multiple collisions, depending on the selections,
11351131
// it may be accepted for some associations and rejected for other
1136-
if (fConfigVariousOptions.fPropMuon) {
1132+
if (static_cast<int>(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) {
11371133
VarManager::FillPropagateMuon<TMuonFillMap>(muon, collision);
11381134
}
1139-
// recalculte pDca and global muon kinematics
1140-
if (static_cast<int>(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) {
1135+
// recalculate pDca / DCA and global muon kinematics
1136+
// kMuonPDca is always taken from MCH (standalone or the MCH matched to a global)
1137+
if (static_cast<int>(muon.trackType()) <= 2) {
11411138
auto muontrack = muon.template matchMCHTrack_as<TMuons>();
1142-
if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) {
1143-
continue;
1144-
}
1145-
auto mfttrack = muon.template matchMFTTrack_as<TMFTTracks>();
11461139
VarManager::FillTrackCollision<TMuonFillMap>(muontrack, collision);
1147-
if constexpr (static_cast<bool>(TMFTFillMap & VarManager::ObjTypes::MFTCov)) {
1148-
auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]);
1149-
VarManager::FillGlobalMuonRefitCov<TMuonFillMap, TMFTFillMap>(muontrack, mfttrack, collision, mfttrackcov);
1140+
if (fConfigVariousOptions.fRefitGlobalMuon) {
1141+
if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) {
1142+
continue;
1143+
}
1144+
auto mfttrack = muon.template matchMFTTrack_as<TMFTTracks>();
1145+
// Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov)
1146+
if constexpr (static_cast<bool>(TMFTFillMap & VarManager::ObjTypes::MFTCov)) {
1147+
auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]);
1148+
VarManager::FillGlobalMuonRefitCov<TMuonFillMap, TMFTFillMap>(muontrack, mfttrack, collision, mfttrackcov);
1149+
} else {
1150+
VarManager::FillGlobalMuonRefit<TMuonFillMap>(muontrack, mfttrack, collision);
1151+
}
11501152
} else {
1151-
VarManager::FillGlobalMuonRefit<TMuonFillMap>(muontrack, mfttrack, collision);
1153+
// Helix DCA of the global track; leaves kMuonPDca from the matched MCH above
1154+
VarManager::FillTrackCollision<TMuonFillMap>(muon, collision);
11521155
}
11531156
} else {
11541157
VarManager::FillTrackCollision<TMuonFillMap>(muon, collision);
@@ -1262,21 +1265,28 @@ struct TableMakerMC {
12621265
}
12631266

12641267
VarManager::FillTrack<TMuonFillMap>(muon);
1265-
if (fConfigVariousOptions.fPropMuon) {
1268+
if (static_cast<int>(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) {
12661269
VarManager::FillPropagateMuon<TMuonFillMap>(muon, collision);
12671270
}
1268-
// recalculte pDca and global muon kinematics
1271+
// recalculate pDca / DCA and global muon kinematics
1272+
// kMuonPDca is always taken from MCH (standalone or the MCH matched to a global)
12691273
int globalClusters = muon.nClusters();
1270-
if (static_cast<int>(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) {
1274+
if (static_cast<int>(muon.trackType()) <= 2) {
12711275
auto muontrack = muon.template matchMCHTrack_as<TMuons>();
1272-
auto mfttrack = muon.template matchMFTTrack_as<TMFTTracks>();
1273-
globalClusters += mfttrack.nClusters();
12741276
VarManager::FillTrackCollision<TMuonFillMap>(muontrack, collision);
1275-
if constexpr (static_cast<bool>(TMFTFillMap & VarManager::ObjTypes::MFTCov)) {
1276-
auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]);
1277-
VarManager::FillGlobalMuonRefitCov<TMuonFillMap, TMFTFillMap>(muontrack, mfttrack, collision, mfttrackcov);
1277+
if (fConfigVariousOptions.fRefitGlobalMuon) {
1278+
auto mfttrack = muon.template matchMFTTrack_as<TMFTTracks>();
1279+
globalClusters += mfttrack.nClusters();
1280+
// Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov)
1281+
if constexpr (static_cast<bool>(TMFTFillMap & VarManager::ObjTypes::MFTCov)) {
1282+
auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]);
1283+
VarManager::FillGlobalMuonRefitCov<TMuonFillMap, TMFTFillMap>(muontrack, mfttrack, collision, mfttrackcov);
1284+
} else {
1285+
VarManager::FillGlobalMuonRefit<TMuonFillMap>(muontrack, mfttrack, collision);
1286+
}
12781287
} else {
1279-
VarManager::FillGlobalMuonRefit<TMuonFillMap>(muontrack, mfttrack, collision);
1288+
// Helix DCA of the global track; leaves kMuonPDca from the matched MCH above
1289+
VarManager::FillTrackCollision<TMuonFillMap>(muon, collision);
12801290
}
12811291
} else {
12821292
VarManager::FillTrackCollision<TMuonFillMap>(muon, collision);
@@ -1323,16 +1333,6 @@ struct TableMakerMC {
13231333
o2::base::Propagator::initFieldFromGRP(fGrpMag);
13241334
VarManager::SetMagneticField(fGrpMag->getNominalL3Field());
13251335
}
1326-
if (fConfigCCDB.fUseRemoteZShift) {
1327-
auto* fZShift = fCCDB->getForTimeStamp<std::vector<float>>(fConfigCCDB.fZShiftPath, bcs.begin().timestamp());
1328-
if (fZShift != nullptr && !fZShift->empty()) {
1329-
VarManager::SetZShift((*fZShift)[0]);
1330-
} else {
1331-
LOG(fatal) << "Could not retrieve Z-shift value from CCDB";
1332-
}
1333-
} else {
1334-
VarManager::SetZShift(fConfigCCDB.fManualZShift.value);
1335-
}
13361336
if (fConfigVariousOptions.fPropMuon) {
13371337
VarManager::SetupMuonMagField();
13381338
}

0 commit comments

Comments
 (0)