diff --git a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/CheckResidConfig.h b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/CheckResidConfig.h index 53dffeed7ad69..2a07eaf87930f 100644 --- a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/CheckResidConfig.h +++ b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/CheckResidConfig.h @@ -26,11 +26,14 @@ struct CheckResidConfig : o2::conf::ConfigurableParamHelper { bool pvcontribOnly = true; bool addPVAsCluster = true; - bool refitPV = true; bool useStableRef = true; bool doIBOB = true; bool doResid = true; + bool refitPV = true; + float refitPVMV = false; + float refitPVIniScale = 100.f; + O2ParamDef(CheckResidConfig, "checkresid"); }; } // namespace o2::checkresid diff --git a/Detectors/GlobalTrackingWorkflow/study/src/CheckResid.cxx b/Detectors/GlobalTrackingWorkflow/study/src/CheckResid.cxx index d665af5747c60..e6584a7055446 100644 --- a/Detectors/GlobalTrackingWorkflow/study/src/CheckResid.cxx +++ b/Detectors/GlobalTrackingWorkflow/study/src/CheckResid.cxx @@ -38,14 +38,17 @@ #include "CommonUtils/TreeStreamRedirector.h" #include "ReconstructionDataFormats/VtxTrackRef.h" #include "DetectorsVertexing/PVertexer.h" - #ifdef WITH_OPENMP #include #endif +// Attention: in case the residuals are checked with geometry different from the one used for initial reconstruction, +// pass a --configKeyValues option for vertex refit as: +// ;pvertexer.useMeanVertexConstraint=false;pvertexer.iniScale2=100;pvertexer.acceptableScale2=10.; +// In any case, it is better to pass ;pvertexer.useMeanVertexConstraint=false; + namespace o2::checkresid { - using namespace o2::framework; using DetID = o2::detectors::DetID; using DataRequest = o2::globaltracking::DataRequest; @@ -83,6 +86,7 @@ class CheckResidSpec : public Task o2::globaltracking::RecoContainer* mRecoData = nullptr; int mNThreads = 1; + bool mMeanVertexUpdated = false; float mITSROFrameLengthMUS = 0.f; o2::dataformats::MeanVertexObject mMeanVtx{}; std::vector> mITSClustersArray; ///< ITS clusters created in run() method from compact clusters @@ -131,6 +135,7 @@ void CheckResidSpec::updateTimeDependentParams(ProcessingContext& pc) // mTPCCorrMapsLoader.extractCCDBInputs(pc); static bool initOnceDone = false; if (!initOnceDone) { // this params need to be queried only once + const auto& params = o2::checkresid::CheckResidConfig::Instance(); initOnceDone = true; // Note: reading of the ITS AlpideParam needed for ITS timing is done by the RecoContainer auto grp = o2::base::GRPGeomHelper::instance().getGRPECS(); @@ -142,9 +147,13 @@ void CheckResidSpec::updateTimeDependentParams(ProcessingContext& pc) } auto geom = o2::its::GeometryTGeo::Instance(); geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G, o2::math_utils::TransformType::T2G)); - o2::conf::ConfigurableParam::updateFromString("pvertexer.useMeanVertexConstraint=false"); + o2::conf::ConfigurableParam::updateFromString("pvertexer.useTimeInChi2=false;"); mVertexer.init(); } + if (mMeanVertexUpdated) { + mMeanVertexUpdated = false; + mVertexer.initMeanVertexConstraint(); + } bool updateMaps = false; /* if (mTPCCorrMapsLoader.isUpdated()) { @@ -200,6 +209,7 @@ void CheckResidSpec::process() } nvGood++; if (params.refitPV) { + LOGP(debug, "Refitting PV#{} of {} tracks", iv, pve.getNContributors()); auto tStartPVF = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); bool res = refitPV(pve, iv); pvFitDuration += std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count() - tStartPVF; @@ -315,6 +325,7 @@ bool CheckResidSpec::processITSTrack(const o2::its::TrackITS& iTrack, const o2:: resTrack.points.clear(); if (!prop->propagateToDCA(pv, trFitOut, bz)) { + LOGP(debug, "Failed to propagateToDCA, {}", trFitOut.asString()); return false; } float cosAlp, sinAlp; @@ -418,7 +429,7 @@ bool CheckResidSpec::refitPV(o2::dataformats::PrimaryVertex& pv, int vid) std::vector tracks; std::vector useTrack; std::vector gidsITS; - int ntr = pv.getNContributors(); + int ntr = pv.getNContributors(), ntrIni = ntr; tracks.reserve(ntr); useTrack.reserve(ntr); gidsITS.reserve(ntr); @@ -447,20 +458,18 @@ bool CheckResidSpec::refitPV(o2::dataformats::PrimaryVertex& pv, int vid) ntr++; } if (ntr < params.minPVContributors || !mVertexer.prepareVertexRefit(tracks, pv)) { + LOGP(warn, "Abandon vertex refit: NcontribNew = {} vs NcontribOld = {}", ntr, ntrIni); return false; } - // readjust vertexZ - const auto& pool = mVertexer.getTracksPool(); - float zUpd = 0; - for (const auto& t : pool) { - zUpd += t.z; - } - if (pool.size()) { - pv.setZ(zUpd / pool.size()); - mVertexer.prepareVertexRefit(tracks, pv); + LOGP(debug, "Original vtx: Nc:{} {}, chi2={}", pv.getNContributors(), pv.asString(), pv.getChi2()); + auto pvSave = pv; + pv = mVertexer.refitVertexFull(useTrack, pv); + LOGP(debug, "Refitted vtx: Nc:{} {}, chi2={}", ntr, pv.asString(), pv.getChi2()); + if (pv.getChi2() < 0.f) { + LOGP(warn, "Failed to refit PV {}", pvSave.asString()); + return false; } - pv = mVertexer.refitVertex(useTrack, pv); - return pv.getChi2() > 0.f; + return true; } bool CheckResidSpec::refitITStrack(o2::track::TrackParCov& track, GTrackID gid) @@ -515,6 +524,7 @@ void CheckResidSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) if (matcher == ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) { LOG(info) << "Imposing new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString(); mMeanVtx = *(const o2::dataformats::MeanVertexObject*)obj; + mMeanVertexUpdated = true; return; } if (matcher == ConcreteDataMatcher("ITS", "CLUSDICT", 0)) { diff --git a/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h b/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h index 9967cbfcd5642..c06c2119b0cd1 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h @@ -113,7 +113,7 @@ class PVertexer bool prepareVertexRefit(const TR& tracks, const o2d::VertexBase& vtxSeed); PVertex refitVertex(const std::vector useTrack, const o2d::VertexBase& vtxSeed); - + PVertex refitVertexFull(const std::vector useTrack, const o2d::VertexBase& vtxSeed); auto getNTZClusters() const { return mNTZClustersIni; } auto getTotTrials() const { return mTotTrials; } auto getMaxTrialsPerCluster() const { return mMaxTrialPerCluster; } @@ -135,6 +135,7 @@ class PVertexer void setPoolDumpDirectory(const std::string& d) { mPoolDumpDirectory = d; } void printInpuTracksStatus(const VertexingInput& input) const; + void initMeanVertexConstraint(); private: static constexpr int DBS_UNDEF = -2, DBS_NOISE = -1, DBS_INCHECK = -10; @@ -152,7 +153,6 @@ class PVertexer FitStatus evalIterations(VertexSeed& vtxSeed, PVertex& vtx) const; TimeEst timeEstimate(const VertexingInput& input) const; float findZSeedHistoPeak() const; - void initMeanVertexConstraint(); void applyConstraint(VertexSeed& vtxSeed) const; bool upscaleSigma(VertexSeed& vtxSeed) const; bool relateTrackToMeanVertex(o2::track::TrackParCov& trc, float vtxErr2); diff --git a/Detectors/Vertexing/src/PVertexer.cxx b/Detectors/Vertexing/src/PVertexer.cxx index 5fea1943ac762..10e504bba0772 100644 --- a/Detectors/Vertexing/src/PVertexer.cxx +++ b/Detectors/Vertexing/src/PVertexer.cxx @@ -1333,6 +1333,37 @@ PVertex PVertexer::refitVertex(const std::vector useTrack, const o2d::Vert return vtxRes; } +//______________________________________________ +PVertex PVertexer::refitVertexFull(const std::vector useTrack, const o2d::VertexBase& vtxSeed) +{ + // Use this method if because of e.g. different alingnment the new vertex is supposed to be shifted from the original one. + // Refit the tracks prepared by the successful prepareVertexRefit, possible skipping those tracks wich have useTrack value false + // (useTrack is ignored if empty). + // The vtxSeed is the originally found vertex, assumed to be the same original PV used for the prepareVertexRefit. + // Refitted PrimaryVertex is returned, negative chi2 means failure of the refit. + // ATTENTION: only the position is refitted, the vertex time and IRMin/IRMax info is dummy. + + if (vtxSeed != mVtxRefitOrig) { + throw std::runtime_error("refitVertex must be preceded by successful prepareVertexRefit"); + } + VertexingInput inp; + inp.scaleSigma2 = mPVParams->iniScale2; + inp.idRange = gsl::span(mRefitTrackIDs); + if (useTrack.size()) { + for (uint32_t i = 0; i < mTracksPool.size(); i++) { + mTracksPool[i].vtxID = useTrack[mTracksPool[i].entry] ? TrackVF::kNoVtx : TrackVF::kDiscarded; + } + } + PVertex vtxRes{}; + vtxRes.VertexBase::operator=(vtxSeed); + if (findVertex(inp, vtxRes)) { + vtxRes.setTimeStamp({0.f, -1.}); // time is not refitter + } else { + vtxRes.setChi2(-1.); + } + return vtxRes; +} + //______________________________________________ void PVertexer::printInpuTracksStatus(const VertexingInput& input) const {