Skip to content

Commit fa60244

Browse files
dsekihatdsekihat
authored
[PWGEM/Dilepton] update createResolutionMap.cxx and studyMCTruth.cxx (#10310)
1 parent a560935 commit fa60244

File tree

2 files changed

+217
-44
lines changed

2 files changed

+217
-44
lines changed

‎PWGEM/Dilepton/Tasks/createResolutionMap.cxx‎

Lines changed: 91 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@
3535

3636
#include "CCDB/BasicCCDBManager.h"
3737
#include "DataFormatsParameters/GRPMagField.h"
38+
#include "DataFormatsCalibration/MeanVertexObject.h"
3839
#include "TGeoGlobalMagField.h"
3940
#include "Field/MagneticField.h"
4041

@@ -66,6 +67,11 @@ struct CreateResolutionMap {
6667
Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
6768
Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
6869
Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
70+
Configurable<std::string> lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"};
71+
Configurable<std::string> mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"};
72+
Configurable<float> d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"};
73+
Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
74+
6975
Configurable<int> cfgEventGeneratorType{"cfgEventGeneratorType", -1, "if positive, select event generator type. i.e. gap or signal"};
7076
Configurable<int> cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"};
7177
Configurable<bool> cfg_require_true_mc_collision_association{"cfg_require_true_mc_collision_association", true, "flag to require true mc collision association"};
@@ -148,7 +154,12 @@ struct CreateResolutionMap {
148154
Service<o2::ccdb::BasicCCDBManager> ccdb;
149155
o2::globaltracking::MatchGlobalFwd mMatching;
150156
int mRunNumber = 0;
151-
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
157+
float d_bz;
158+
// o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
159+
o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
160+
o2::dataformats::VertexBase mVtx;
161+
const o2::dataformats::MeanVertexObject* mMeanVtx = nullptr;
162+
o2::base::MatLayerCylSet* lut = nullptr;
152163

153164
void init(o2::framework::InitContext&)
154165
{
@@ -166,6 +177,9 @@ struct CreateResolutionMap {
166177
ccdb->setFatalWhenNull(false);
167178
ccdbApi.init(ccdburl);
168179

180+
mRunNumber = 0;
181+
d_bz = 0;
182+
169183
const AxisSpec axis_cent{ConfCentBins, "centrality (%)"};
170184
const AxisSpec axis_pt_gen{ConfPtGenBins, "p_{T,l}^{gen} (GeV/c)"};
171185
const AxisSpec axis_eta_cb_gen{ConfEtaCBGenBins, "#eta_{l}^{gen}"};
@@ -199,12 +213,67 @@ struct CreateResolutionMap {
199213
return;
200214
}
201215

216+
// load matLUT for this timestamp
217+
if (!lut) {
218+
LOG(info) << "Loading material look-up table for timestamp: " << bc.timestamp();
219+
lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp<o2::base::MatLayerCylSet>(lutPath, bc.timestamp()));
220+
} else {
221+
LOG(info) << "Material look-up table already in place. Not reloading.";
222+
}
223+
224+
// In case override, don't proceed, please - no CCDB access required
225+
if (d_bz_input > -990) {
226+
d_bz = d_bz_input;
227+
o2::parameters::GRPMagField grpmag;
228+
if (std::fabs(d_bz) > 1e-5) {
229+
grpmag.setL3Current(30000.f / (d_bz / 5.0f));
230+
}
231+
o2::base::Propagator::initFieldFromGRP(&grpmag);
232+
o2::base::Propagator::Instance()->setMatLUT(lut);
233+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
234+
mRunNumber = bc.runNumber();
235+
236+
if (!o2::base::GeometryManager::isGeometryLoaded()) {
237+
ccdb->get<TGeoManager>(geoPath);
238+
}
239+
o2::mch::TrackExtrap::setField();
240+
return;
241+
}
242+
243+
auto run3grp_timestamp = bc.timestamp();
244+
o2::parameters::GRPObject* grpo = 0x0;
245+
o2::parameters::GRPMagField* grpmag = 0x0;
246+
if (!skipGRPOquery) {
247+
grpo = ccdb->getForTimeStamp<o2::parameters::GRPObject>(grpPath, run3grp_timestamp);
248+
}
249+
if (grpo) {
250+
o2::base::Propagator::initFieldFromGRP(grpo);
251+
o2::base::Propagator::Instance()->setMatLUT(lut);
252+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
253+
// Fetch magnetic field from ccdb for current collision
254+
d_bz = grpo->getNominalL3Field();
255+
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
256+
} else {
257+
grpmag = ccdb->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, run3grp_timestamp);
258+
if (!grpmag) {
259+
LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp;
260+
}
261+
o2::base::Propagator::initFieldFromGRP(grpmag);
262+
o2::base::Propagator::Instance()->setMatLUT(lut);
263+
mMeanVtx = ccdb->getForTimeStamp<o2::dataformats::MeanVertexObject>(mVtxPath, bc.timestamp());
264+
265+
// Fetch magnetic field from ccdb for current collision
266+
d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f);
267+
LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG";
268+
}
202269
mRunNumber = bc.runNumber();
203-
std::map<string, string> metadata;
204-
auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, mRunNumber);
205-
auto ts = soreor.first;
206-
auto grpmag = ccdbApi.retrieveFromTFileAny<o2::parameters::GRPMagField>(grpmagPath, metadata, ts);
207-
o2::base::Propagator::initFieldFromGRP(grpmag);
270+
271+
// std::map<string, string> metadata;
272+
// auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, mRunNumber);
273+
// auto ts = soreor.first;
274+
// auto grpmag = ccdbApi.retrieveFromTFileAny<o2::parameters::GRPMagField>(grpmagPath, metadata, ts);
275+
// o2::base::Propagator::initFieldFromGRP(grpmag);
276+
208277
if (!o2::base::GeometryManager::isGeometryLoaded()) {
209278
ccdb->get<TGeoManager>(geoPath);
210279
}
@@ -337,16 +406,19 @@ struct CreateResolutionMap {
337406
return false;
338407
}
339408

340-
gpu::gpustd::array<float, 2> dcaInfo;
409+
o2::dataformats::DCA mDcaInfoCov;
410+
mDcaInfoCov.set(999, 999, 999, 999, 999);
341411
auto track_par_cov_recalc = getTrackParCov(track);
342412
track_par_cov_recalc.setPID(o2::track::PID::Electron);
343-
std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
344-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
345-
getPxPyPz(track_par_cov_recalc, pVec_recalc);
413+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
414+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
415+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
416+
float dcaXY = mDcaInfoCov.getY();
417+
float dcaZ = mDcaInfoCov.getZ();
346418

347419
// LOGF(info, "collision.globalIndex() = %d, track.collisionId() = %d, track.pt() = %.16f, track_par_cov_recalc.getPt() = %.16f", collision.globalIndex(), track.collisionId(), track.pt(), track_par_cov_recalc.getPt());
348420

349-
if (std::fabs(dcaInfo[0]) > electroncuts.cfg_max_dcaxy || std::fabs(dcaInfo[1]) > electroncuts.cfg_max_dcaz) {
421+
if (std::fabs(dcaXY) > electroncuts.cfg_max_dcaxy || std::fabs(dcaZ) > electroncuts.cfg_max_dcaz) {
350422
return false;
351423
}
352424

@@ -412,7 +484,7 @@ struct CreateResolutionMap {
412484
void fillMuon(TCollision const& collision, TMuon const& muon, const float centrality)
413485
{
414486
auto mcparticle = muon.template mcParticle_as<aod::McParticles>();
415-
if (abs(mcparticle.pdgCode()) != 13 || !(mcparticle.isPhysicalPrimary() || mcparticle.producedByGenerator())) {
487+
if (std::abs(mcparticle.pdgCode()) != 13 || !(mcparticle.isPhysicalPrimary() || mcparticle.producedByGenerator())) {
416488
return;
417489
}
418490
if (cfg_require_true_mc_collision_association && mcparticle.mcCollisionId() != collision.mcCollisionId()) {
@@ -527,12 +599,15 @@ struct CreateResolutionMap {
527599
return;
528600
}
529601

530-
gpu::gpustd::array<float, 2> dcaInfo;
602+
o2::dataformats::DCA mDcaInfoCov;
603+
mDcaInfoCov.set(999, 999, 999, 999, 999);
531604
auto track_par_cov_recalc = getTrackParCov(track);
532605
track_par_cov_recalc.setPID(o2::track::PID::Electron);
533-
// std::array<float, 3> pVec_recalc = {0, 0, 0}; // px, py, pz
534-
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo);
535-
// getPxPyPz(track_par_cov_recalc, pVec_recalc);
606+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
607+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
608+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, track_par_cov_recalc, 2.f, matCorr, &mDcaInfoCov);
609+
// float dcaXY = mDcaInfoCov.getY();
610+
// float dcaZ = mDcaInfoCov.getZ();
536611

537612
float pt = track_par_cov_recalc.getPt();
538613
float eta = track_par_cov_recalc.getEta();

0 commit comments

Comments
 (0)