Skip to content

Commit 5443f10

Browse files
dsekihatdsekihat
authored
[PWGEM/Dilepton] add collision cov. matrix and matlut (#10242)
1 parent c72059c commit 5443f10

File tree

1 file changed

+96
-31
lines changed

1 file changed

+96
-31
lines changed

‎PWGEM/Dilepton/TableProducer/skimmerPrimaryElectron.cxx‎

Lines changed: 96 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include "DetectorsBase/GeometryManager.h"
2626
#include "DataFormatsParameters/GRPObject.h"
2727
#include "DataFormatsParameters/GRPMagField.h"
28+
#include "DataFormatsCalibration/MeanVertexObject.h"
2829
#include "CCDB/BasicCCDBManager.h"
2930
#include "Common/Core/trackUtilities.h"
3031
#include "CommonConstants/PhysicsConstants.h"
@@ -62,6 +63,8 @@ struct skimmerPrimaryElectron {
6263
Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
6364
Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
6465
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
66+
Configurable<std::string> lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"};
67+
Configurable<std::string> mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"};
6568
Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
6669

6770
// Operation and minimisation criteria
@@ -97,7 +100,11 @@ struct skimmerPrimaryElectron {
97100
int mRunNumber;
98101
float d_bz;
99102
Service<o2::ccdb::BasicCCDBManager> ccdb;
100-
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
103+
// o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
104+
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
105+
o2::dataformats::VertexBase mVtx;
106+
const o2::dataformats::MeanVertexObject* mMeanVtx = nullptr;
107+
o2::base::MatLayerCylSet* lut = nullptr;
101108

102109
void init(InitContext&)
103110
{
@@ -154,6 +161,14 @@ struct skimmerPrimaryElectron {
154161
return;
155162
}
156163

164+
// load matLUT for this timestamp
165+
if (!lut) {
166+
LOG(info) << "Loading material look-up table for timestamp: " << bc.timestamp();
167+
lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp<o2::base::MatLayerCylSet>(lutPath, bc.timestamp()));
168+
} else {
169+
LOG(info) << "Material look-up table already in place. Not reloading.";
170+
}
171+
157172
// In case override, don't proceed, please - no CCDB access required
158173
if (d_bz_input > -990) {
159174
d_bz = d_bz_input;
@@ -162,17 +177,22 @@ struct skimmerPrimaryElectron {
162177
grpmag.setL3Current(30000.f / (d_bz / 5.0f));
163178
}
164179
o2::base::Propagator::initFieldFromGRP(&grpmag);
180+
o2::base::Propagator::Instance()->setMatLUT(lut);
181+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
165182
mRunNumber = bc.runNumber();
166183
return;
167184
}
168185

169186
auto run3grp_timestamp = bc.timestamp();
170187
o2::parameters::GRPObject* grpo = 0x0;
171188
o2::parameters::GRPMagField* grpmag = 0x0;
172-
if (!skipGRPOquery)
189+
if (!skipGRPOquery) {
173190
grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(grpPath, run3grp_timestamp);
191+
}
174192
if (grpo) {
175193
o2::base::Propagator::initFieldFromGRP(grpo);
194+
o2::base::Propagator::Instance()->setMatLUT(lut);
195+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
176196
// Fetch magnetic field from ccdb for current collision
177197
d_bz = grpo->getNominalL3Field();
178198
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
@@ -182,6 +202,9 @@ struct skimmerPrimaryElectron {
182202
LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp;
183203
}
184204
o2::base::Propagator::initFieldFromGRP(grpmag);
205+
o2::base::Propagator::Instance()->setMatLUT(lut);
206+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
207+
185208
// Fetch magnetic field from ccdb for current collision
186209
d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f);
187210
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
@@ -240,14 +263,17 @@ struct skimmerPrimaryElectron {
240263
return false;
241264
}
242265

243-
gpu::gpustd::array<float, 2> dcaInfo;
266+
o2::dataformats::DCA mDcaInfoCov;
267+
mDcaInfoCov.set(999, 999, 999, 999, 999);
244268
auto track_par_cov_recalc = getTrackParCov(track);
245269
track_par_cov_recalc.setPID(o2::track::PID::Electron);
246-
// std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
247-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
248-
// getPxPyPz(track_par_cov_recalc, pVec_recalc);
249-
float dcaXY = dcaInfo[0];
250-
float dcaZ = dcaInfo[1];
270+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
271+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
272+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
273+
float dcaXY = mDcaInfoCov.getY();
274+
float dcaZ = mDcaInfoCov.getZ();
275+
276+
// LOGF(info, "track_par_cov_recalc.getSigmaY2() = %.16f, mDcaInfoCov.getSigmaY2() = %.16f, track_par_cov_recalc.getSigmaZ2() = %.16f, mDcaInfoCov.getSigmaZ2() = %.16f, track_par_cov_recalc.getSigmaZY() = %.16f, mDcaInfoCov.getSigmaYZ() = %.16f", track_par_cov_recalc.getSigmaY2(), mDcaInfoCov.getSigmaY2(), track_par_cov_recalc.getSigmaZ2(), mDcaInfoCov.getSigmaZ2(), track_par_cov_recalc.getSigmaZY(), mDcaInfoCov.getSigmaYZ());
251277

252278
if (std::fabs(dcaXY) > dca_xy_max || std::fabs(dcaZ) > dca_z_max) {
253279
return false;
@@ -312,14 +338,15 @@ struct skimmerPrimaryElectron {
312338
void fillTrackTable(TCollision const& collision, TTrack const& track)
313339
{
314340
if (std::find(stored_trackIds.begin(), stored_trackIds.end(), std::pair<int, int>{collision.globalIndex(), track.globalIndex()}) == stored_trackIds.end()) {
315-
gpu::gpustd::array<float, 2> dcaInfo;
341+
o2::dataformats::DCA mDcaInfoCov;
342+
mDcaInfoCov.set(999, 999, 999, 999, 999);
316343
auto track_par_cov_recalc = getTrackParCov(track);
317344
track_par_cov_recalc.setPID(o2::track::PID::Electron);
318-
// std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
319-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
320-
// getPxPyPz(track_par_cov_recalc, pVec_recalc);
321-
float dcaXY = dcaInfo[0];
322-
float dcaZ = dcaInfo[1];
345+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
346+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
347+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
348+
float dcaXY = mDcaInfoCov.getY();
349+
float dcaZ = mDcaInfoCov.getZ();
323350

324351
float pt_recalc = track_par_cov_recalc.getPt();
325352
float eta_recalc = track_par_cov_recalc.getEta();
@@ -620,6 +647,8 @@ struct prefilterPrimaryElectron {
620647
Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
621648
Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
622649
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
650+
Configurable<std::string> lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"};
651+
Configurable<std::string> mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"};
623652
Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
624653

625654
// Operation and minimisation criteria
@@ -649,7 +678,11 @@ struct prefilterPrimaryElectron {
649678
int mRunNumber;
650679
float d_bz;
651680
Service<o2::ccdb::BasicCCDBManager> ccdb;
652-
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
681+
// o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
682+
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
683+
o2::dataformats::VertexBase mVtx;
684+
const o2::dataformats::MeanVertexObject* mMeanVtx = nullptr;
685+
o2::base::MatLayerCylSet* lut = nullptr;
653686

654687
void init(InitContext&)
655688
{
@@ -684,6 +717,14 @@ struct prefilterPrimaryElectron {
684717
return;
685718
}
686719

720+
// load matLUT for this timestamp
721+
if (!lut) {
722+
LOG(info) << "Loading material look-up table for timestamp: " << bc.timestamp();
723+
lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp<o2::base::MatLayerCylSet>(lutPath, bc.timestamp()));
724+
} else {
725+
LOG(info) << "Material look-up table already in place. Not reloading.";
726+
}
727+
687728
// In case override, don't proceed, please - no CCDB access required
688729
if (d_bz_input > -990) {
689730
d_bz = d_bz_input;
@@ -692,6 +733,8 @@ struct prefilterPrimaryElectron {
692733
grpmag.setL3Current(30000.f / (d_bz / 5.0f));
693734
}
694735
o2::base::Propagator::initFieldFromGRP(&grpmag);
736+
o2::base::Propagator::Instance()->setMatLUT(lut);
737+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
695738
mRunNumber = bc.runNumber();
696739
return;
697740
}
@@ -703,6 +746,8 @@ struct prefilterPrimaryElectron {
703746
grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(grpPath, run3grp_timestamp);
704747
if (grpo) {
705748
o2::base::Propagator::initFieldFromGRP(grpo);
749+
o2::base::Propagator::Instance()->setMatLUT(lut);
750+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
706751
// Fetch magnetic field from ccdb for current collision
707752
d_bz = grpo->getNominalL3Field();
708753
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
@@ -712,6 +757,8 @@ struct prefilterPrimaryElectron {
712757
LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp;
713758
}
714759
o2::base::Propagator::initFieldFromGRP(grpmag);
760+
o2::base::Propagator::Instance()->setMatLUT(lut);
761+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
715762
// Fetch magnetic field from ccdb for current collision
716763
d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f);
717764
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
@@ -759,14 +806,17 @@ struct prefilterPrimaryElectron {
759806
return false;
760807
}
761808

762-
gpu::gpustd::array<float, 2> dcaInfo;
809+
o2::dataformats::DCA mDcaInfoCov;
810+
mDcaInfoCov.set(999, 999, 999, 999, 999);
763811
auto track_par_cov_recalc = getTrackParCov(track);
764812
track_par_cov_recalc.setPID(o2::track::PID::Electron);
765-
// std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
766-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
767-
// getPxPyPz(track_par_cov_recalc, pVec_recalc);
813+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
814+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
815+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
816+
float dcaXY = mDcaInfoCov.getY();
817+
float dcaZ = mDcaInfoCov.getZ();
768818

769-
if (std::fabs(dcaInfo[0]) > max_dcaxy || std::fabs(dcaInfo[1]) > max_dcaz) {
819+
if (std::fabs(dcaXY) > max_dcaxy || std::fabs(dcaZ) > max_dcaz) {
770820
return false;
771821
}
772822

@@ -781,12 +831,15 @@ struct prefilterPrimaryElectron {
781831
bool reconstructPC(TCollision const& collision, TTrack1 const& ele, TTrack2 const& pos)
782832
{
783833
float mee = 0, phiv = 0;
784-
gpu::gpustd::array<float, 2> dcaInfo;
834+
o2::dataformats::DCA mDcaInfoCov;
835+
mDcaInfoCov.set(999, 999, 999, 999, 999);
785836
std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
837+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
838+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
786839

787840
if constexpr (loose_track_sign > 0) { // positive track is loose track
788841
auto track_par_cov_recalc = getTrackParCov(pos);
789-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
842+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
790843
getPxPyPz(track_par_cov_recalc, pVec_recalc);
791844

792845
ROOT::Math::PtEtaPhiMVector v1(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron);
@@ -796,7 +849,7 @@ struct prefilterPrimaryElectron {
796849
phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pVec_recalc[0], pVec_recalc[1], pVec_recalc[2], ele.px(), ele.py(), ele.pz(), pos.sign(), ele.sign(), d_bz);
797850
} else {
798851
auto track_par_cov_recalc = getTrackParCov(ele);
799-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
852+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
800853
getPxPyPz(track_par_cov_recalc, pVec_recalc);
801854

802855
ROOT::Math::PtEtaPhiMVector v1(track_par_cov_recalc.getPt(), track_par_cov_recalc.getEta(), track_par_cov_recalc.getPhi(), o2::constants::physics::MassElectron);
@@ -861,11 +914,14 @@ struct prefilterPrimaryElectron {
861914
if (!checkTrack(collision, ele)) {
862915
continue;
863916
}
864-
gpu::gpustd::array<float, 2> dcaInfo;
917+
o2::dataformats::DCA mDcaInfoCov;
918+
mDcaInfoCov.set(999, 999, 999, 999, 999);
865919
std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
866920
auto track_par_cov_recalc = getTrackParCov(ele);
867921
track_par_cov_recalc.setPID(o2::track::PID::Electron);
868-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
922+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
923+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
924+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
869925
getPxPyPz(track_par_cov_recalc, pVec_recalc);
870926

871927
for (auto& empos : positrons_per_coll) {
@@ -899,11 +955,14 @@ struct prefilterPrimaryElectron {
899955
if (!checkTrack(collision, pos)) { // track cut is applied to loose sample
900956
continue;
901957
}
902-
gpu::gpustd::array<float, 2> dcaInfo;
958+
o2::dataformats::DCA mDcaInfoCov;
959+
mDcaInfoCov.set(999, 999, 999, 999, 999);
903960
std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
904961
auto track_par_cov_recalc = getTrackParCov(pos);
905962
track_par_cov_recalc.setPID(o2::track::PID::Electron);
906-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
963+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
964+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
965+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
907966
getPxPyPz(track_par_cov_recalc, pVec_recalc);
908967
for (auto& emele : electrons_per_coll) {
909968
if (emele.trackId() == pos.globalIndex()) {
@@ -935,11 +994,14 @@ struct prefilterPrimaryElectron {
935994
if (!checkTrack(collision, pos)) { // track cut is applied to loose sample
936995
continue;
937996
}
938-
gpu::gpustd::array<float, 2> dcaInfo;
997+
o2::dataformats::DCA mDcaInfoCov;
998+
mDcaInfoCov.set(999, 999, 999, 999, 999);
939999
std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
9401000
auto track_par_cov_recalc = getTrackParCov(pos);
9411001
track_par_cov_recalc.setPID(o2::track::PID::Electron);
942-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
1002+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
1003+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
1004+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
9431005
getPxPyPz(track_par_cov_recalc, pVec_recalc);
9441006
for (auto& empos : positrons_per_coll) {
9451007
if (empos.trackId() == pos.globalIndex()) {
@@ -959,11 +1021,14 @@ struct prefilterPrimaryElectron {
9591021
if (!checkTrack(collision, ele)) {
9601022
continue;
9611023
}
962-
gpu::gpustd::array<float, 2> dcaInfo;
1024+
o2::dataformats::DCA mDcaInfoCov;
1025+
mDcaInfoCov.set(999, 999, 999, 999, 999);
9631026
std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
9641027
auto track_par_cov_recalc = getTrackParCov(ele);
9651028
track_par_cov_recalc.setPID(o2::track::PID::Electron);
966-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
1029+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
1030+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
1031+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
9671032
getPxPyPz(track_par_cov_recalc, pVec_recalc);
9681033

9691034
for (auto& emele : electrons_per_coll) {

0 commit comments

Comments
 (0)