Skip to content

[PWGDQ] Add function to propagate muon tracks through absorber to PV#3762

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
lucamicheletti93 merged 13 commits into from
Nov 20, 2023
Merged

[PWGDQ] Add function to propagate muon tracks through absorber to PV #3762

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
3dd984f
Add function to propagate muon tracks through absorber to PV
mcoquet642 Nov 2, 2023
ad0825b
Format
mcoquet642 Nov 2, 2023
f22bc42
Fix
mcoquet642 Nov 3, 2023
8cc8dfb
Fix
mcoquet642 Nov 7, 2023
6c2ae57
Fix typo
mcoquet642 Nov 7, 2023
e6f2130
Fixing libraries
mcoquet642 Nov 9, 2023
14c905d
Clang format
mcoquet642 Nov 9, 2023
f987b36
Missing array definition
mcoquet642 Nov 9, 2023
1d597e2
Retrieving geometry from ccdb
mcoquet642 Nov 10, 2023
ff2991e
Clang format
mcoquet642 Nov 10, 2023
b2fab9e
Fixing retrieving of ccdb object
mcoquet642 Nov 15, 2023
d7d43ef
Calng format
mcoquet642 Nov 15, 2023
aa80bad
Fixing covariance name
mcoquet642 Nov 17, 2023
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
2 changes: 1 addition & 1 deletion PWGDQ/Core/CMakeLists.txt
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -21,7 +21,7 @@ o2physics_add_library(PWGDQCore
AnalysisCompositeCut.cxx
MCProng.cxx
MCSignal.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore KFParticle::KFParticle)
PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2::GlobalTracking O2Physics::AnalysisCore KFParticle::KFParticle)

o2physics_target_root_dictionary(PWGDQCore
HEADERS AnalysisCut.h
Expand Down
1 change: 1 addition & 0 deletions PWGDQ/Core/VarManager.cxx
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -33,6 +33,7 @@ o2::vertexing::DCAFitterN<2> VarManager::fgFitterTwoProngBarrel;
o2::vertexing::DCAFitterN<3> VarManager::fgFitterThreeProngBarrel;
o2::vertexing::FwdDCAFitterN<2> VarManager::fgFitterTwoProngFwd;
o2::vertexing::FwdDCAFitterN<3> VarManager::fgFitterThreeProngFwd;
o2::globaltracking::MatchGlobalFwd VarManager::mMatching;
std::map<VarManager::CalibObjects, TObject*> VarManager::fgCalibs;
bool VarManager::fgRunTPCPostCalibration[4] = {false, false, false, false};

Expand Down
80 changes: 80 additions & 0 deletions PWGDQ/Core/VarManager.h
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -36,6 +36,7 @@
#include "Math/GenVector/Boost.h"

#include "Framework/DataTypes.h"
// #include "MCHTracking/TrackExtrap.h"
#include "ReconstructionDataFormats/Track.h"
#include "ReconstructionDataFormats/Vertex.h"
#include "DCAFitter/DCAFitterN.h"
Expand All@@ -46,6 +47,7 @@
#include "Math/SMatrix.h"
#include "ReconstructionDataFormats/TrackFwd.h"
#include "DCAFitter/FwdDCAFitterN.h"
#include "GlobalTracking/MatchGlobalFwd.h"
#include "CommonConstants/PhysicsConstants.h"

#include "KFParticle.h"
Expand DownExpand Up@@ -194,9 +196,13 @@ class VarManager : public TObject
kNEventWiseVariables,

// Basic track/muon/pair wise variables
kX,
kY,
kZ,
kPt,
kInvPt,
kEta,
kTgl,
kPhi,
kP,
kPx,
Expand DownExpand Up@@ -309,9 +315,19 @@ class VarManager : public TObject
kMuonChi2MatchMCHMFT,
kMuonMatchScoreMCHMFT,
kMuonCXX,
kMuonCXY,
kMuonCYY,
kMuonCPhiX,
kMuonCPhiY,
kMuonCPhiPhi,
kMuonCTglX,
kMuonCTglY,
kMuonCTglPhi,
kMuonCTglTgl,
kMuonC1Pt2X,
kMuonC1Pt2Y,
kMuonC1Pt2Phi,
kMuonC1Pt2Tgl,
kMuonC1Pt21Pt2,
kNMuonTrackVariables,
kMuonTrackType,
Expand DownExpand Up@@ -558,6 +574,8 @@ class VarManager : public TObject
return (1.0 / harm) * TMath::ATan(qnya / qnxa);
};

template <uint32_t fillMap, typename T, typename C>
static void FillPropagateMuon(const T& muon, const C& collision, float* values = nullptr);
template <uint32_t fillMap, typename T>
static void FillEvent(T const& event, float* values = nullptr);
template <uint32_t fillMap, typename T>
Expand DownExpand Up@@ -649,6 +667,7 @@ class VarManager : public TObject
static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel;
static o2::vertexing::FwdDCAFitterN<2> fgFitterTwoProngFwd;
static o2::vertexing::FwdDCAFitterN<3> fgFitterThreeProngFwd;
static o2::globaltracking::MatchGlobalFwd mMatching;

static std::map<CalibObjects, TObject*> fgCalibs; // map of calibration histograms
static bool fgRunTPCPostCalibration[4]; // 0-electron, 1-pion, 2-kaon, 3-proton
Expand DownExpand Up@@ -743,6 +762,53 @@ KFPVertex VarManager::createKFPVertexFromCollision(const T& collision)
return kfpVertex;
}

template <uint32_t fillMap, typename T, typename C>
void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* values)
{
if (!values) {
values = fgValues;
}
if constexpr ((fillMap & MuonCov) > 0) {
o2::mch::TrackExtrap::setField();
double chi2 = muon.chi2();
SMatrix5 tpars(muon.x(), muon.y(), muon.phi(), muon.tgl(), muon.signed1Pt());
std::vector<double> v1{muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(),
muon.cPhiPhi(), muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(),
muon.c1PtX(), muon.c1PtY(), muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2()};
SMatrix55 tcovs(v1.begin(), v1.end());
o2::track::TrackParCovFwd fwdtrack{muon.z(), tpars, tcovs, chi2};
o2::dataformats::GlobalFwdTrack track{fwdtrack};
auto mchTrack = mMatching.FwdtoMCH(track);
o2::mch::TrackExtrap::extrapToVertex(mchTrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY());
auto propmuon = mMatching.MCHtoFwd(mchTrack);

values[kPt] = propmuon.getPt();
values[kX] = propmuon.getX();
values[kY] = propmuon.getY();
values[kZ] = propmuon.getZ();
values[kEta] = propmuon.getEta();
values[kTgl] = propmuon.getTgl();
values[kPhi] = propmuon.getPhi();

SMatrix55 cov = propmuon.getCovariances();
values[kMuonCXX] = cov(0, 0);
values[kMuonCXY] = cov(1, 0);
values[kMuonCYY] = cov(1, 1);
values[kMuonCPhiX] = cov(2, 0);
values[kMuonCPhiY] = cov(2, 1);
values[kMuonCPhiPhi] = cov(2, 2);
values[kMuonCTglX] = cov(3, 0);
values[kMuonCTglY] = cov(3, 1);
values[kMuonCTglPhi] = cov(3, 2);
values[kMuonCTglTgl] = cov(3, 3);
values[kMuonC1Pt2X] = cov(4, 0);
values[kMuonC1Pt2Y] = cov(4, 1);
values[kMuonC1Pt2Phi] = cov(4, 2);
values[kMuonC1Pt2Tgl] = cov(4, 3);
values[kMuonC1Pt21Pt2] = cov(4, 4);
}
}

template <uint32_t fillMap, typename T>
void VarManager::FillEvent(T const& event, float* values)
{
Expand DownExpand Up@@ -1268,10 +1334,24 @@ void VarManager::FillTrack(T const& track, float* values)
}
// Quantities based on the muon covariance table
if constexpr ((fillMap & ReducedMuonCov) > 0 || (fillMap & MuonCov) > 0) {
values[kX] = track.x();
values[kY] = track.y();
values[kZ] = track.z();
values[kTgl] = track.tgl();
values[kMuonCXX] = track.cXX();
values[kMuonCXY] = track.cXY();
values[kMuonCYY] = track.cYY();
values[kMuonCPhiX] = track.cPhiX();
values[kMuonCPhiY] = track.cPhiY();
values[kMuonCPhiPhi] = track.cPhiPhi();
values[kMuonCTglX] = track.cTglX();
values[kMuonCTglY] = track.cTglY();
values[kMuonCTglPhi] = track.cTglPhi();
values[kMuonCTglTgl] = track.cTglTgl();
values[kMuonC1Pt2X] = track.c1PtX();
values[kMuonC1Pt2Y] = track.c1PtY();
values[kMuonC1Pt2Phi] = track.c1PtPhi();
values[kMuonC1Pt2Tgl] = track.c1PtTgl();
values[kMuonC1Pt21Pt2] = track.c1Pt21Pt2();
}

Expand Down
74 changes: 53 additions & 21 deletions PWGDQ/TableProducer/tableMaker.cxx
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -49,6 +49,11 @@
#include "MathUtils/Primitive2D.h"
#include "DataFormatsGlobalTracking/RecoContainer.h"
#include "Common/DataModel/CollisionAssociationTables.h"
#include "DataFormatsParameters/GRPMagField.h"
#include "Field/MagneticField.h"
#include "TGeoGlobalMagField.h"
#include "DetectorsBase/Propagator.h"
#include "DetectorsBase/GeometryManager.h"

using std::cout;
using std::endl;
Expand DownExpand Up@@ -163,8 +168,14 @@ struct TableMaker {
Configurable<bool> fConfigSaveElectronSample{"cfgSaveElectronSample", false, "If true, only save electron sample"};
Configurable<bool> fConfigDummyRunlist{"cfgDummyRunlist", false, "If true, use dummy runlist"};
Configurable<int> fConfigInitRunNumber{"cfgInitRunNumber", 543215, "Initial run number used in run by run checks"};
Configurable<bool> fPropMuon{"cfgPropMuon", false, "Propgate muon tracks through absorber"};
Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};

Service<o2::ccdb::BasicCCDBManager> fCCDB;

o2::parameters::GRPMagField* grpmag = nullptr;

AnalysisCompositeCut* fEventCut; //! Event selection cut
std::vector<AnalysisCompositeCut> fTrackCuts; //! Barrel track cuts
std::vector<AnalysisCompositeCut> fMuonCuts; //! Muon track cuts
Expand All@@ -174,6 +185,7 @@ struct TableMaker {
Preslice<aod::TrackAssoc> trackIndicesPerCollision = aod::track_association::collisionId;
Preslice<aod::FwdTrackAssoc> fwdtrackIndicesPerCollision = aod::track_association::collisionId;

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

Expand All@@ -185,6 +197,14 @@ struct TableMaker {
void init(o2::framework::InitContext& context)
{
DefineCuts();
ccdb->setURL(fConfigCcdbUrl);
ccdb->setCaching(true);
ccdb->setLocalObjectValidityChecking();
if (fPropMuon) {
if (!o2::base::GeometryManager::isGeometryLoaded()) {
ccdb->get<TGeoManager>(geoPath);
}
}

VarManager::SetDefaultVarNames();
fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars);
Expand DownExpand Up@@ -319,17 +339,23 @@ struct TableMaker {
void fullSkimming(TEvent const& collision, aod::BCsWithTimestamps const&, TTracks const& tracksBarrel, TMuons const& tracksMuon, TAmbiTracks const& ambiTracksMid, TAmbiMuons const& ambiTracksFwd, TMFTTracks const& mftTracks = nullptr)
{
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
if (fConfigComputeTPCpostCalib && fCurrentRun != bc.runNumber()) {
auto calibList = fCCDB->getForTimeStamp<TList>(fConfigCcdbPathTPC.value, bc.timestamp());
VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron"));
VarManager::SetCalibrationObject(VarManager::kTPCElectronSigma, calibList->FindObject("sigma_map_electron"));
VarManager::SetCalibrationObject(VarManager::kTPCPionMean, calibList->FindObject("mean_map_pion"));
VarManager::SetCalibrationObject(VarManager::kTPCPionSigma, calibList->FindObject("sigma_map_pion"));
VarManager::SetCalibrationObject(VarManager::kTPCProtonMean, calibList->FindObject("mean_map_proton"));
VarManager::SetCalibrationObject(VarManager::kTPCProtonSigma, calibList->FindObject("sigma_map_proton"));
if (fConfigComputeTPCpostCalibKaon) {
VarManager::SetCalibrationObject(VarManager::kTPCKaonMean, calibList->FindObject("mean_map_kaon"));
VarManager::SetCalibrationObject(VarManager::kTPCKaonSigma, calibList->FindObject("sigma_map_kaon"));
if (fCurrentRun != bc.runNumber()) {
if (fConfigComputeTPCpostCalib) {
auto calibList = fCCDB->getForTimeStamp<TList>(fConfigCcdbPathTPC.value, bc.timestamp());
VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron"));
VarManager::SetCalibrationObject(VarManager::kTPCElectronSigma, calibList->FindObject("sigma_map_electron"));
VarManager::SetCalibrationObject(VarManager::kTPCPionMean, calibList->FindObject("mean_map_pion"));
VarManager::SetCalibrationObject(VarManager::kTPCPionSigma, calibList->FindObject("sigma_map_pion"));
VarManager::SetCalibrationObject(VarManager::kTPCProtonMean, calibList->FindObject("mean_map_proton"));
VarManager::SetCalibrationObject(VarManager::kTPCProtonSigma, calibList->FindObject("sigma_map_proton"));
if (fConfigComputeTPCpostCalibKaon) {
VarManager::SetCalibrationObject(VarManager::kTPCKaonMean, calibList->FindObject("mean_map_kaon"));
VarManager::SetCalibrationObject(VarManager::kTPCKaonSigma, calibList->FindObject("sigma_map_kaon"));
}
}
grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, bc.timestamp());
if (grpmag != nullptr) {
o2::base::Propagator::initFieldFromGRP(grpmag);
}
fCurrentRun = bc.runNumber();
}
Expand DownExpand Up@@ -603,6 +629,9 @@ struct TableMaker {
trackTempFilterMap = uint8_t(0);

VarManager::FillTrack<TMuonFillMap>(muon);
if (fPropMuon) {
VarManager::FillPropagateMuon<TMuonFillMap>(muon, collision);
}
if (fDoDetailedQA) {
fHistMan->FillHistClass("Muons_BeforeCuts", VarManager::fgValues);
if (fIsAmbiguous && isAmbiguous == 1) {
Expand DownExpand Up@@ -663,18 +692,18 @@ struct TableMaker {
}
}

muonBasic(event.lastIndex(), trackFilteringTag, muon.pt(), muon.eta(), muon.phi(), muon.sign(), isAmbiguous);
muonBasic(event.lastIndex(), trackFilteringTag, VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], muon.sign(), isAmbiguous);
muonExtra(muon.nClusters(), muon.pDca(), muon.rAtAbsorberEnd(),
muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(),
muon.matchScoreMCHMFT(), newMatchIndex.find(muon.index())->second, newMFTMatchIndex.find(muon.index())->second, muon.mchBitMap(), muon.midBitMap(),
muon.midBoards(), muon.trackType(), muon.fwdDcaX(), muon.fwdDcaY(),
muon.trackTime(), muon.trackTimeRes());
muonInfo(muon.collisionId(), collision.posX(), collision.posY(), collision.posZ());
if constexpr (static_cast<bool>(TMuonFillMap & VarManager::ObjTypes::MuonCov)) {
muonCov(muon.x(), muon.y(), muon.z(), muon.phi(), muon.tgl(), muon.signed1Pt(),
muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), muon.cPhiPhi(),
muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), muon.c1PtX(), muon.c1PtY(),
muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2());
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],
VarManager::fgValues[VarManager::kMuonCXX], VarManager::fgValues[VarManager::kMuonCXY], VarManager::fgValues[VarManager::kMuonCYY], VarManager::fgValues[VarManager::kMuonCPhiX], VarManager::fgValues[VarManager::kMuonCPhiY], VarManager::fgValues[VarManager::kMuonCPhiPhi],
VarManager::fgValues[VarManager::kMuonCTglX], VarManager::fgValues[VarManager::kMuonCTglY], VarManager::fgValues[VarManager::kMuonCTglPhi], VarManager::fgValues[VarManager::kMuonCTglTgl], VarManager::fgValues[VarManager::kMuonC1Pt2X], VarManager::fgValues[VarManager::kMuonC1Pt2Y],
VarManager::fgValues[VarManager::kMuonC1Pt2Phi], VarManager::fgValues[VarManager::kMuonC1Pt2Tgl], VarManager::fgValues[VarManager::kMuonC1Pt21Pt2]);
}
}
} // end if constexpr (TMuonFillMap)
Expand DownExpand Up@@ -924,6 +953,9 @@ struct TableMaker {
trackTempFilterMap = uint8_t(0);

VarManager::FillTrack<TMuonFillMap>(muon);
if (fPropMuon) {
VarManager::FillPropagateMuon<TMuonFillMap>(muon, collision);
}
if (fDoDetailedQA) {
fHistMan->FillHistClass("Muons_BeforeCuts", VarManager::fgValues);
if (fIsAmbiguous && isAmbiguous == 1) {
Expand DownExpand Up@@ -970,18 +1002,18 @@ struct TableMaker {
}
}

muonBasic(event.lastIndex(), trackFilteringTag, muon.pt(), muon.eta(), muon.phi(), muon.sign(), isAmbiguous);
muonBasic(event.lastIndex(), trackFilteringTag, VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], muon.sign(), isAmbiguous);
muonExtra(muon.nClusters(), muon.pDca(), muon.rAtAbsorberEnd(),
muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(),
muon.matchScoreMCHMFT(), newMatchIndex.find(muon.index())->second, -1, muon.mchBitMap(), muon.midBitMap(),
muon.midBoards(), muon.trackType(), muon.fwdDcaX(), muon.fwdDcaY(),
muon.trackTime(), muon.trackTimeRes());
muonInfo(muon.collisionId(), collision.posX(), collision.posY(), collision.posZ());
if constexpr (static_cast<bool>(TMuonFillMap & VarManager::ObjTypes::MuonCov)) {
muonCov(muon.x(), muon.y(), muon.z(), muon.phi(), muon.tgl(), muon.signed1Pt(),
muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), muon.cPhiPhi(),
muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), muon.c1PtX(), muon.c1PtY(),
muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2());
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],
VarManager::fgValues[VarManager::kMuonCXX], VarManager::fgValues[VarManager::kMuonCXY], VarManager::fgValues[VarManager::kMuonCYY], VarManager::fgValues[VarManager::kMuonCPhiX], VarManager::fgValues[VarManager::kMuonCPhiY], VarManager::fgValues[VarManager::kMuonCPhiPhi],
VarManager::fgValues[VarManager::kMuonCTglX], VarManager::fgValues[VarManager::kMuonCTglY], VarManager::fgValues[VarManager::kMuonCTglPhi], VarManager::fgValues[VarManager::kMuonCTglTgl], VarManager::fgValues[VarManager::kMuonC1Pt2X], VarManager::fgValues[VarManager::kMuonC1Pt2Y],
VarManager::fgValues[VarManager::kMuonC1Pt2Phi], VarManager::fgValues[VarManager::kMuonC1Pt2Tgl], VarManager::fgValues[VarManager::kMuonC1Pt21Pt2]);
}
}
} // end if constexpr (TMuonFillMap)
Expand Down