Skip to content

Commit 89c76d3

Browse files
mcoquet642mcoquet642
authored
[PWGDQ] Add function to propagate muon tracks through absorber to PV (#3762)
* Add function to propagate muon tracks through absorber to PV * Format * Fix * Fix * Fix typo * Fixing libraries * Clang format * Missing array definition * Retrieving geometry from ccdb * Clang format * Fixing retrieving of ccdb object * Calng format * Fixing covariance name
1 parent 303ed96 commit 89c76d3

File tree

4 files changed

+135
-22
lines changed

4 files changed

+135
-22
lines changed

‎PWGDQ/Core/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ o2physics_add_library(PWGDQCore
2121
AnalysisCompositeCut.cxx
2222
MCProng.cxx
2323
MCSignal.cxx
24-
PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore KFParticle::KFParticle)
24+
PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2::GlobalTracking O2Physics::AnalysisCore KFParticle::KFParticle)
2525

2626
o2physics_target_root_dictionary(PWGDQCore
2727
HEADERS AnalysisCut.h

‎PWGDQ/Core/VarManager.cxx‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ o2::vertexing::DCAFitterN<2> VarManager::fgFitterTwoProngBarrel;
3333
o2::vertexing::DCAFitterN<3> VarManager::fgFitterThreeProngBarrel;
3434
o2::vertexing::FwdDCAFitterN<2> VarManager::fgFitterTwoProngFwd;
3535
o2::vertexing::FwdDCAFitterN<3> VarManager::fgFitterThreeProngFwd;
36+
o2::globaltracking::MatchGlobalFwd VarManager::mMatching;
3637
std::map<VarManager::CalibObjects, TObject*> VarManager::fgCalibs;
3738
bool VarManager::fgRunTPCPostCalibration[4] = {false, false, false, false};
3839

‎PWGDQ/Core/VarManager.h‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
#include "Math/GenVector/Boost.h"
3737

3838
#include "Framework/DataTypes.h"
39+
// #include "MCHTracking/TrackExtrap.h"
3940
#include "ReconstructionDataFormats/Track.h"
4041
#include "ReconstructionDataFormats/Vertex.h"
4142
#include "DCAFitter/DCAFitterN.h"
@@ -46,6 +47,7 @@
4647
#include "Math/SMatrix.h"
4748
#include "ReconstructionDataFormats/TrackFwd.h"
4849
#include "DCAFitter/FwdDCAFitterN.h"
50+
#include "GlobalTracking/MatchGlobalFwd.h"
4951
#include "CommonConstants/PhysicsConstants.h"
5052

5153
#include "KFParticle.h"
@@ -194,10 +196,14 @@ class VarManager : public TObject
194196
kNEventWiseVariables,
195197

196198
// Basic track/muon/pair wise variables
199+
kX,
200+
kY,
201+
kZ,
197202
kPt,
198203
kSignedPt,
199204
kInvPt,
200205
kEta,
206+
kTgl,
201207
kPhi,
202208
kP,
203209
kPx,
@@ -310,9 +316,19 @@ class VarManager : public TObject
310316
kMuonChi2MatchMCHMFT,
311317
kMuonMatchScoreMCHMFT,
312318
kMuonCXX,
319+
kMuonCXY,
313320
kMuonCYY,
321+
kMuonCPhiX,
322+
kMuonCPhiY,
314323
kMuonCPhiPhi,
324+
kMuonCTglX,
325+
kMuonCTglY,
326+
kMuonCTglPhi,
315327
kMuonCTglTgl,
328+
kMuonC1Pt2X,
329+
kMuonC1Pt2Y,
330+
kMuonC1Pt2Phi,
331+
kMuonC1Pt2Tgl,
316332
kMuonC1Pt21Pt2,
317333
kNMuonTrackVariables,
318334
kMuonTrackType,
@@ -559,6 +575,8 @@ class VarManager : public TObject
559575
return (1.0 / harm) * TMath::ATan(qnya / qnxa);
560576
};
561577

578+
template <uint32_t fillMap, typename T, typename C>
579+
static void FillPropagateMuon(const T& muon, const C& collision, float* values = nullptr);
562580
template <uint32_t fillMap, typename T>
563581
static void FillEvent(T const& event, float* values = nullptr);
564582
template <uint32_t fillMap, typename T>
@@ -650,6 +668,7 @@ class VarManager : public TObject
650668
static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel;
651669
static o2::vertexing::FwdDCAFitterN<2> fgFitterTwoProngFwd;
652670
static o2::vertexing::FwdDCAFitterN<3> fgFitterThreeProngFwd;
671+
static o2::globaltracking::MatchGlobalFwd mMatching;
653672

654673
static std::map<CalibObjects, TObject*> fgCalibs; // map of calibration histograms
655674
static bool fgRunTPCPostCalibration[4]; // 0-electron, 1-pion, 2-kaon, 3-proton
@@ -744,6 +763,53 @@ KFPVertex VarManager::createKFPVertexFromCollision(const T& collision)
744763
return kfpVertex;
745764
}
746765

766+
template <uint32_t fillMap, typename T, typename C>
767+
void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* values)
768+
{
769+
if (!values) {
770+
values = fgValues;
771+
}
772+
if constexpr ((fillMap & MuonCov) > 0) {
773+
o2::mch::TrackExtrap::setField();
774+
double chi2 = muon.chi2();
775+
SMatrix5 tpars(muon.x(), muon.y(), muon.phi(), muon.tgl(), muon.signed1Pt());
776+
std::vector<double> v1{muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(),
777+
muon.cPhiPhi(), muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(),
778+
muon.c1PtX(), muon.c1PtY(), muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2()};
779+
SMatrix55 tcovs(v1.begin(), v1.end());
780+
o2::track::TrackParCovFwd fwdtrack{muon.z(), tpars, tcovs, chi2};
781+
o2::dataformats::GlobalFwdTrack track{fwdtrack};
782+
auto mchTrack = mMatching.FwdtoMCH(track);
783+
o2::mch::TrackExtrap::extrapToVertex(mchTrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY());
784+
auto propmuon = mMatching.MCHtoFwd(mchTrack);
785+
786+
values[kPt] = propmuon.getPt();
787+
values[kX] = propmuon.getX();
788+
values[kY] = propmuon.getY();
789+
values[kZ] = propmuon.getZ();
790+
values[kEta] = propmuon.getEta();
791+
values[kTgl] = propmuon.getTgl();
792+
values[kPhi] = propmuon.getPhi();
793+
794+
SMatrix55 cov = propmuon.getCovariances();
795+
values[kMuonCXX] = cov(0, 0);
796+
values[kMuonCXY] = cov(1, 0);
797+
values[kMuonCYY] = cov(1, 1);
798+
values[kMuonCPhiX] = cov(2, 0);
799+
values[kMuonCPhiY] = cov(2, 1);
800+
values[kMuonCPhiPhi] = cov(2, 2);
801+
values[kMuonCTglX] = cov(3, 0);
802+
values[kMuonCTglY] = cov(3, 1);
803+
values[kMuonCTglPhi] = cov(3, 2);
804+
values[kMuonCTglTgl] = cov(3, 3);
805+
values[kMuonC1Pt2X] = cov(4, 0);
806+
values[kMuonC1Pt2Y] = cov(4, 1);
807+
values[kMuonC1Pt2Phi] = cov(4, 2);
808+
values[kMuonC1Pt2Tgl] = cov(4, 3);
809+
values[kMuonC1Pt21Pt2] = cov(4, 4);
810+
}
811+
}
812+
747813
template <uint32_t fillMap, typename T>
748814
void VarManager::FillEvent(T const& event, float* values)
749815
{
@@ -1270,10 +1336,24 @@ void VarManager::FillTrack(T const& track, float* values)
12701336
}
12711337
// Quantities based on the muon covariance table
12721338
if constexpr ((fillMap & ReducedMuonCov) > 0 || (fillMap & MuonCov) > 0) {
1339+
values[kX] = track.x();
1340+
values[kY] = track.y();
1341+
values[kZ] = track.z();
1342+
values[kTgl] = track.tgl();
12731343
values[kMuonCXX] = track.cXX();
1344+
values[kMuonCXY] = track.cXY();
12741345
values[kMuonCYY] = track.cYY();
1346+
values[kMuonCPhiX] = track.cPhiX();
1347+
values[kMuonCPhiY] = track.cPhiY();
12751348
values[kMuonCPhiPhi] = track.cPhiPhi();
1349+
values[kMuonCTglX] = track.cTglX();
1350+
values[kMuonCTglY] = track.cTglY();
1351+
values[kMuonCTglPhi] = track.cTglPhi();
12761352
values[kMuonCTglTgl] = track.cTglTgl();
1353+
values[kMuonC1Pt2X] = track.c1PtX();
1354+
values[kMuonC1Pt2Y] = track.c1PtY();
1355+
values[kMuonC1Pt2Phi] = track.c1PtPhi();
1356+
values[kMuonC1Pt2Tgl] = track.c1PtTgl();
12771357
values[kMuonC1Pt21Pt2] = track.c1Pt21Pt2();
12781358
}
12791359

‎PWGDQ/TableProducer/tableMaker.cxx‎

Lines changed: 53 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,11 @@
4949
#include "MathUtils/Primitive2D.h"
5050
#include "DataFormatsGlobalTracking/RecoContainer.h"
5151
#include "Common/DataModel/CollisionAssociationTables.h"
52+
#include "DataFormatsParameters/GRPMagField.h"
53+
#include "Field/MagneticField.h"
54+
#include "TGeoGlobalMagField.h"
55+
#include "DetectorsBase/Propagator.h"
56+
#include "DetectorsBase/GeometryManager.h"
5257

5358
using std::cout;
5459
using std::endl;
@@ -163,8 +168,14 @@ struct TableMaker {
163168
Configurable<bool> fConfigSaveElectronSample{"cfgSaveElectronSample", false, "If true, only save electron sample"};
164169
Configurable<bool> fConfigDummyRunlist{"cfgDummyRunlist", false, "If true, use dummy runlist"};
165170
Configurable<int> fConfigInitRunNumber{"cfgInitRunNumber", 543215, "Initial run number used in run by run checks"};
171+
Configurable<bool> fPropMuon{"cfgPropMuon", false, "Propgate muon tracks through absorber"};
172+
Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
173+
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
174+
166175
Service<o2::ccdb::BasicCCDBManager> fCCDB;
167176

177+
o2::parameters::GRPMagField* grpmag = nullptr;
178+
168179
AnalysisCompositeCut* fEventCut; //! Event selection cut
169180
std::vector<AnalysisCompositeCut> fTrackCuts; //! Barrel track cuts
170181
std::vector<AnalysisCompositeCut> fMuonCuts; //! Muon track cuts
@@ -174,6 +185,7 @@ struct TableMaker {
174185
Preslice<aod::TrackAssoc> trackIndicesPerCollision = aod::track_association::collisionId;
175186
Preslice<aod::FwdTrackAssoc> fwdtrackIndicesPerCollision = aod::track_association::collisionId;
176187

188+
Service<o2::ccdb::BasicCCDBManager> ccdb;
177189
bool fDoDetailedQA = false; // Bool to set detailed QA true, if QA is set true
178190
int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc.
179191

@@ -185,6 +197,14 @@ struct TableMaker {
185197
void init(o2::framework::InitContext& context)
186198
{
187199
DefineCuts();
200+
ccdb->setURL(fConfigCcdbUrl);
201+
ccdb->setCaching(true);
202+
ccdb->setLocalObjectValidityChecking();
203+
if (fPropMuon) {
204+
if (!o2::base::GeometryManager::isGeometryLoaded()) {
205+
ccdb->get<TGeoManager>(geoPath);
206+
}
207+
}
188208

189209
VarManager::SetDefaultVarNames();
190210
fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars);
@@ -319,17 +339,23 @@ struct TableMaker {
319339
void fullSkimming(TEvent const& collision, aod::BCsWithTimestamps const&, TTracks const& tracksBarrel, TMuons const& tracksMuon, TAmbiTracks const& ambiTracksMid, TAmbiMuons const& ambiTracksFwd, TMFTTracks const& mftTracks = nullptr)
320340
{
321341
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
322-
if (fConfigComputeTPCpostCalib && fCurrentRun != bc.runNumber()) {
323-
auto calibList = fCCDB->getForTimeStamp<TList>(fConfigCcdbPathTPC.value, bc.timestamp());
324-
VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron"));
325-
VarManager::SetCalibrationObject(VarManager::kTPCElectronSigma, calibList->FindObject("sigma_map_electron"));
326-
VarManager::SetCalibrationObject(VarManager::kTPCPionMean, calibList->FindObject("mean_map_pion"));
327-
VarManager::SetCalibrationObject(VarManager::kTPCPionSigma, calibList->FindObject("sigma_map_pion"));
328-
VarManager::SetCalibrationObject(VarManager::kTPCProtonMean, calibList->FindObject("mean_map_proton"));
329-
VarManager::SetCalibrationObject(VarManager::kTPCProtonSigma, calibList->FindObject("sigma_map_proton"));
330-
if (fConfigComputeTPCpostCalibKaon) {
331-
VarManager::SetCalibrationObject(VarManager::kTPCKaonMean, calibList->FindObject("mean_map_kaon"));
332-
VarManager::SetCalibrationObject(VarManager::kTPCKaonSigma, calibList->FindObject("sigma_map_kaon"));
342+
if (fCurrentRun != bc.runNumber()) {
343+
if (fConfigComputeTPCpostCalib) {
344+
auto calibList = fCCDB->getForTimeStamp<TList>(fConfigCcdbPathTPC.value, bc.timestamp());
345+
VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron"));
346+
VarManager::SetCalibrationObject(VarManager::kTPCElectronSigma, calibList->FindObject("sigma_map_electron"));
347+
VarManager::SetCalibrationObject(VarManager::kTPCPionMean, calibList->FindObject("mean_map_pion"));
348+
VarManager::SetCalibrationObject(VarManager::kTPCPionSigma, calibList->FindObject("sigma_map_pion"));
349+
VarManager::SetCalibrationObject(VarManager::kTPCProtonMean, calibList->FindObject("mean_map_proton"));
350+
VarManager::SetCalibrationObject(VarManager::kTPCProtonSigma, calibList->FindObject("sigma_map_proton"));
351+
if (fConfigComputeTPCpostCalibKaon) {
352+
VarManager::SetCalibrationObject(VarManager::kTPCKaonMean, calibList->FindObject("mean_map_kaon"));
353+
VarManager::SetCalibrationObject(VarManager::kTPCKaonSigma, calibList->FindObject("sigma_map_kaon"));
354+
}
355+
}
356+
grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, bc.timestamp());
357+
if (grpmag != nullptr) {
358+
o2::base::Propagator::initFieldFromGRP(grpmag);
333359
}
334360
fCurrentRun = bc.runNumber();
335361
}
@@ -603,6 +629,9 @@ struct TableMaker {
603629
trackTempFilterMap = uint8_t(0);
604630

605631
VarManager::FillTrack<TMuonFillMap>(muon);
632+
if (fPropMuon) {
633+
VarManager::FillPropagateMuon<TMuonFillMap>(muon, collision);
634+
}
606635
if (fDoDetailedQA) {
607636
fHistMan->FillHistClass("Muons_BeforeCuts", VarManager::fgValues);
608637
if (fIsAmbiguous && isAmbiguous == 1) {
@@ -663,18 +692,18 @@ struct TableMaker {
663692
}
664693
}
665694

666-
muonBasic(event.lastIndex(), trackFilteringTag, muon.pt(), muon.eta(), muon.phi(), muon.sign(), isAmbiguous);
695+
muonBasic(event.lastIndex(), trackFilteringTag, VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], muon.sign(), isAmbiguous);
667696
muonExtra(muon.nClusters(), muon.pDca(), muon.rAtAbsorberEnd(),
668697
muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(),
669698
muon.matchScoreMCHMFT(), newMatchIndex.find(muon.index())->second, newMFTMatchIndex.find(muon.index())->second, muon.mchBitMap(), muon.midBitMap(),
670699
muon.midBoards(), muon.trackType(), muon.fwdDcaX(), muon.fwdDcaY(),
671700
muon.trackTime(), muon.trackTimeRes());
672701
muonInfo(muon.collisionId(), collision.posX(), collision.posY(), collision.posZ());
673702
if constexpr (static_cast<bool>(TMuonFillMap & VarManager::ObjTypes::MuonCov)) {
674-
muonCov(muon.x(), muon.y(), muon.z(), muon.phi(), muon.tgl(), muon.signed1Pt(),
675-
muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), muon.cPhiPhi(),
676-
muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), muon.c1PtX(), muon.c1PtY(),
677-
muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2());
703+
muonCov(VarManager::fgValues[VarManager::kX], VarManager::fgValues[VarManager::kY], VarManager::fgValues[VarManager::kZ], VarManager::fgValues[VarManager::kPhi], VarManager::fgValues[VarManager::kTgl], muon.sign() / VarManager::fgValues[VarManager::kPt],
704+
VarManager::fgValues[VarManager::kMuonCXX], VarManager::fgValues[VarManager::kMuonCXY], VarManager::fgValues[VarManager::kMuonCYY], VarManager::fgValues[VarManager::kMuonCPhiX], VarManager::fgValues[VarManager::kMuonCPhiY], VarManager::fgValues[VarManager::kMuonCPhiPhi],
705+
VarManager::fgValues[VarManager::kMuonCTglX], VarManager::fgValues[VarManager::kMuonCTglY], VarManager::fgValues[VarManager::kMuonCTglPhi], VarManager::fgValues[VarManager::kMuonCTglTgl], VarManager::fgValues[VarManager::kMuonC1Pt2X], VarManager::fgValues[VarManager::kMuonC1Pt2Y],
706+
VarManager::fgValues[VarManager::kMuonC1Pt2Phi], VarManager::fgValues[VarManager::kMuonC1Pt2Tgl], VarManager::fgValues[VarManager::kMuonC1Pt21Pt2]);
678707
}
679708
}
680709
} // end if constexpr (TMuonFillMap)
@@ -924,6 +953,9 @@ struct TableMaker {
924953
trackTempFilterMap = uint8_t(0);
925954

926955
VarManager::FillTrack<TMuonFillMap>(muon);
956+
if (fPropMuon) {
957+
VarManager::FillPropagateMuon<TMuonFillMap>(muon, collision);
958+
}
927959
if (fDoDetailedQA) {
928960
fHistMan->FillHistClass("Muons_BeforeCuts", VarManager::fgValues);
929961
if (fIsAmbiguous && isAmbiguous == 1) {
@@ -970,18 +1002,18 @@ struct TableMaker {
9701002
}
9711003
}
9721004

973-
muonBasic(event.lastIndex(), trackFilteringTag, muon.pt(), muon.eta(), muon.phi(), muon.sign(), isAmbiguous);
1005+
muonBasic(event.lastIndex(), trackFilteringTag, VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], muon.sign(), isAmbiguous);
9741006
muonExtra(muon.nClusters(), muon.pDca(), muon.rAtAbsorberEnd(),
9751007
muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(),
9761008
muon.matchScoreMCHMFT(), newMatchIndex.find(muon.index())->second, -1, muon.mchBitMap(), muon.midBitMap(),
9771009
muon.midBoards(), muon.trackType(), muon.fwdDcaX(), muon.fwdDcaY(),
9781010
muon.trackTime(), muon.trackTimeRes());
9791011
muonInfo(muon.collisionId(), collision.posX(), collision.posY(), collision.posZ());
9801012
if constexpr (static_cast<bool>(TMuonFillMap & VarManager::ObjTypes::MuonCov)) {
981-
muonCov(muon.x(), muon.y(), muon.z(), muon.phi(), muon.tgl(), muon.signed1Pt(),
982-
muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), muon.cPhiPhi(),
983-
muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), muon.c1PtX(), muon.c1PtY(),
984-
muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2());
1013+
muonCov(VarManager::fgValues[VarManager::kX], VarManager::fgValues[VarManager::kY], VarManager::fgValues[VarManager::kZ], VarManager::fgValues[VarManager::kPhi], VarManager::fgValues[VarManager::kTgl], muon.sign() / VarManager::fgValues[VarManager::kPt],
1014+
VarManager::fgValues[VarManager::kMuonCXX], VarManager::fgValues[VarManager::kMuonCXY], VarManager::fgValues[VarManager::kMuonCYY], VarManager::fgValues[VarManager::kMuonCPhiX], VarManager::fgValues[VarManager::kMuonCPhiY], VarManager::fgValues[VarManager::kMuonCPhiPhi],
1015+
VarManager::fgValues[VarManager::kMuonCTglX], VarManager::fgValues[VarManager::kMuonCTglY], VarManager::fgValues[VarManager::kMuonCTglPhi], VarManager::fgValues[VarManager::kMuonCTglTgl], VarManager::fgValues[VarManager::kMuonC1Pt2X], VarManager::fgValues[VarManager::kMuonC1Pt2Y],
1016+
VarManager::fgValues[VarManager::kMuonC1Pt2Phi], VarManager::fgValues[VarManager::kMuonC1Pt2Tgl], VarManager::fgValues[VarManager::kMuonC1Pt21Pt2]);
9851017
}
9861018
}
9871019
} // end if constexpr (TMuonFillMap)

0 commit comments

Comments
 (0)