diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h index bb3ebef84032c..953094b5155b7 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h @@ -23,12 +23,14 @@ #include "DetectorsBase/GRPGeomHelper.h" #include "TPCCalibration/VDriftHelper.h" #include "DataFormatsITSMFT/TopologyDictionary.h" +#include "Steer/MCKinematicsReader.h" using namespace o2::framework; namespace o2::globaltracking { struct DataRequest; +struct RecoContainer; } // namespace o2::globaltracking namespace o2 @@ -40,7 +42,7 @@ class TPCInterpolationDPL : public Task public: TPCInterpolationDPL(std::shared_ptr dr, o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap, std::shared_ptr gr, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput, bool extDetResid) : mDataRequest(dr), mSources(src), mSourcesMap(srcMap), mGGCCDBRequest(gr), mUseMC(useMC), mProcessITSTPConly(processITSTPConly), mSendTrackData(sendTrackData), mDebugOutput(debugOutput), mExtDetResid(extDetResid) {} - ~TPCInterpolationDPL() override = default; + ~TPCInterpolationDPL() override; void init(InitContext& ic) final; void run(ProcessingContext& pc) final; void endOfStream(EndOfStreamContext& ec) final; @@ -48,6 +50,7 @@ class TPCInterpolationDPL : public Task private: void updateTimeDependentParams(ProcessingContext& pc); + void fillMCTruth(const o2::globaltracking::RecoContainer& recoData); o2::tpc::TrackInterpolation mInterpolation; ///< track interpolation engine std::shared_ptr mDataRequest; ///< steers the input std::shared_ptr mGGCCDBRequest; @@ -56,6 +59,8 @@ class TPCInterpolationDPL : public Task o2::dataformats::GlobalTrackID::mask_t mSources{}; ///< which input sources are configured o2::dataformats::GlobalTrackID::mask_t mSourcesMap{}; ///< possible subset of mSources specifically for map creation bool mUseMC{false}; ///< MC flag + std::unique_ptr mMCReader; ///< MC kinematics and track references (MC only) + std::vector mTrackDataMC; ///< MC truth aligned with the TrackData output (MC only) bool mProcessITSTPConly{false}; ///< should also tracks without outer point (ITS-TPC only) be processed? bool mProcessSeeds{false}; ///< process not only most complete track, but also its shorter parts bool mDebugOutput{false}; ///< add more information to the output (track points of ITS, TRD and TOF) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h index 99f20e390a09a..e17e8b49509b4 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualAggregatorSpec.h @@ -41,7 +41,7 @@ namespace calibration class ResidualAggregatorDevice : public o2::framework::Task { public: - ResidualAggregatorDevice(std::shared_ptr req, bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData, std::shared_ptr dataRequest) : mCCDBRequest(req), mTrackInput(trackInput), mCTPInput(ctpInput), mWriteUnbinnedResiduals(writeUnbinnedResiduals), mWriteBinnedResiduals(writeBinnedResiduals), mWriteTrackData(writeTrackData), mDataRequest(dataRequest) {} + ResidualAggregatorDevice(std::shared_ptr req, bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData, bool mcInput, std::shared_ptr dataRequest) : mCCDBRequest(req), mTrackInput(trackInput), mCTPInput(ctpInput), mWriteUnbinnedResiduals(writeUnbinnedResiduals), mWriteBinnedResiduals(writeBinnedResiduals), mWriteTrackData(writeTrackData), mMCInput(mcInput), mDataRequest(dataRequest) {} void init(o2::framework::InitContext& ic) final { @@ -97,6 +97,7 @@ class ResidualAggregatorDevice : public o2::framework::Task mAggregator->setWriteBinnedResiduals(mWriteBinnedResiduals); mAggregator->setWriteUnbinnedResiduals(mWriteUnbinnedResiduals); mAggregator->setWriteTrackData(mWriteTrackData); + mAggregator->setWriteTrackDataMC(mMCInput); mAggregator->setCompression(ic.options().get("compression")); } @@ -141,6 +142,14 @@ class ResidualAggregatorDevice : public o2::framework::Task trkData.emplace(pc.inputs().get>("trkData")); trkDataPtr = &trkData.value(); } + // MC truth of the track data (optional, MC only) + const gsl::span* trkDataMCPtr = nullptr; + using trkDataMCType = std::decay_t>(""))>; + std::optional trkDataMC; + if (mMCInput) { + trkDataMC.emplace(pc.inputs().get>("trkDataMC")); + trkDataMCPtr = &trkDataMC.value(); + } // CTP lumi input (optional) const o2::ctp::LumiInfo* lumi = nullptr; using lumiDataType = std::decay_t(""))>; @@ -152,7 +161,7 @@ class ResidualAggregatorDevice : public o2::framework::Task o2::base::TFIDInfoHelper::fillTFIDInfo(pc, mAggregator->getCurrentTFInfo()); LOG(detail) << "Processing TF " << mAggregator->getCurrentTFInfo().tfCounter << " with " << trkData->size() << " tracks and " << residualsData.size() << " unbinned residuals associated to them"; - mAggregator->process(residualsData, residualsDataDet, trackRefs, trkDataPtr, lumi); + mAggregator->process(residualsData, residualsDataDet, trackRefs, trkDataPtr, trkDataMCPtr, lumi); std::chrono::duration runDuration = std::chrono::high_resolution_clock::now() - runStartTime; LOGP(debug, "Duration for run method: {} ms. From this taken for time dependent param update: {} ms", std::chrono::duration_cast(runDuration).count(), @@ -205,6 +214,7 @@ class ResidualAggregatorDevice : public o2::framework::Task bool mWriteBinnedResiduals{false}; ///< flag, whether to write binned residuals to output file bool mWriteUnbinnedResiduals{false}; ///< flag, whether to write unbinned residuals to output file bool mWriteTrackData{false}; ///< flag, whether to write track data to output file + bool mMCInput{false}; ///< flag whether to expect the MC truth of the track data as input bool mRunStopRequested{false}; ///< flag in case the run was stopped bool mInitDone{false}; ///< flag whether initialization was done for current run }; @@ -214,7 +224,7 @@ class ResidualAggregatorDevice : public o2::framework::Task namespace framework { -DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData) +DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, bool writeUnbinnedResiduals, bool writeBinnedResiduals, bool writeTrackData, bool mcInput = false) { std::shared_ptr dataRequest = std::make_shared(); if (ctpInput) { @@ -227,6 +237,9 @@ DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, b inputs.emplace_back("trackRefs", "GLO", "TRKREFS"); if (trackInput) { inputs.emplace_back("trkData", "GLO", "TRKDATA"); + if (mcInput) { + inputs.emplace_back("trkDataMC", "GLO", "TRKDATAMC"); + } } auto ccdbRequest = std::make_shared(true, // orbitResetTime true, // GRPECS=true @@ -240,7 +253,7 @@ DataProcessorSpec getTPCResidualAggregatorSpec(bool trackInput, bool ctpInput, b "residual-aggregator", inputs, Outputs{}, - AlgorithmSpec{adaptFromTask(ccdbRequest, trackInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData, dataRequest)}, + AlgorithmSpec{adaptFromTask(ccdbRequest, trackInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData, trackInput && mcInput, dataRequest)}, Options{ {"sec-per-slot", VariantType::UInt32, 600u, {"number of seconds per calibration time slot (put 0 for infinite slot length)"}}, {"updateInterval", VariantType::UInt32, 6'000u, {"update interval in number of TFs (only used in case slot length is infinite)"}}, diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx index 2af349be4fd37..b7e514bf76516 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx @@ -13,6 +13,8 @@ #include #include +#include +#include #include "DataFormatsITS/TrackITS.h" #include "ITSBase/GeometryTGeo.h" @@ -33,6 +35,9 @@ #include "Framework/ConfigParamRegistry.h" #include "Framework/ControlService.h" #include "Framework/DeviceSpec.h" +#include "Steer/MCKinematicsReader.h" +#include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/O2DatabasePDG.h" using namespace o2::framework; using namespace o2::globaltracking; @@ -44,6 +49,8 @@ namespace o2 namespace tpc { +TPCInterpolationDPL::~TPCInterpolationDPL() = default; + void TPCInterpolationDPL::init(InitContext& ic) { //-------- init geometry and field --------// @@ -59,6 +66,16 @@ void TPCInterpolationDPL::init(InitContext& ic) int lane = ic.services().get().inputTimesliceId; int maxLanes = ic.services().get().maxInputTimeslices; mInterpolation.setLane(lane, maxLanes); + if (mUseMC) { + if (!mSendTrackData) { + LOG(warning) << "MC truth is stored aligned with the track data, but send-track-data is not set: no MC truth will be sent"; + } + mMCReader = std::make_unique(); + auto mcContext = ic.options().get("mc-collision-context"); + if (!mMCReader->initFromDigitContext(mcContext)) { + LOG(fatal) << "Could not initialize the MC kinematics reader from " << mcContext; + } + } } void TPCInterpolationDPL::updateTimeDependentParams(ProcessingContext& pc) @@ -149,9 +166,255 @@ void TPCInterpolationDPL::run(ProcessingContext& pc) if (mDebugOutput) { pc.outputs().snapshot(Output{"GLO", "TRKDATAEXT", 0}, mInterpolation.getTrackDataExtended()); } + if (mUseMC && mSendTrackData) { + fillMCTruth(recoData); + pc.outputs().snapshot(Output{"GLO", "TRKDATAMC", 0}, mTrackDataMC); + } mInterpolation.reset(); } +void TPCInterpolationDPL::fillMCTruth(const RecoContainer& recoData) +{ + // MC truth for every stored TrackData: labels of the ITS-TPC part of the seed and of its ITS and TPC parts, the truth at + // the ITS outer parameters (the ITS track reference nearest to TrackData::par, propagated to its x with the material + // correction of the workflow and the mass of the true particle), the truth at the TPC entrance (the first TPC track + // reference in time, in the sector frame), for tracks with TRD residuals the truth at the TRD entrance and the true + // positions at the x of the TRD tracklets of the track (ideal tracklets), and the origin of the particle of the ITS-TPC + // part (mother, production vertex and process, other stored tracks with the same mother) + const auto& trkData = mInterpolation.getReferenceTracks(); + mTrackDataMC.clear(); + mTrackDataMC.resize(trkData.size()); + struct Lookup { + o2::MCCompLabel lbl{}; + uint32_t idx{0}; + uint8_t kind{0}; // 0: ITS outer, 1: TPC entrance, 2: TRD, 3: origin + }; + std::vector lookups; + lookups.reserve(2 * trkData.size()); + for (size_t i = 0; i < trkData.size(); ++i) { + auto& mc = mTrackDataMC[i]; + auto gidSet = recoData.getSingleDetectorRefs(trkData[i].gid); + auto gidITS = gidSet[GTrackID::ITS].isIndexSet() ? gidSet[GTrackID::ITS] : gidSet[GTrackID::ITSAB]; + if (gidSet[GTrackID::ITSTPC].isIndexSet()) { + mc.label = recoData.getTrackMCLabel(gidSet[GTrackID::ITSTPC]); + } + if (gidITS.isIndexSet()) { + mc.labelITS = recoData.getTrackMCLabel(gidITS); + } + if (gidSet[GTrackID::TPC].isIndexSet()) { + mc.labelTPC = recoData.getTrackMCLabel(gidSet[GTrackID::TPC]); + } + if (mc.labelITS.isValid() && mc.labelTPC.isValid() && mc.labelITS.getTrackEventSourceID() != mc.labelTPC.getTrackEventSourceID()) { + mc.flags |= TrackDataMC::FakeITSTPC; + } + const auto& lblITS = mc.labelITS.isValid() ? mc.labelITS : mc.label; + const auto& lblTPC = mc.labelTPC.isValid() ? mc.labelTPC : mc.label; + if (lblITS.isValid()) { + lookups.push_back({lblITS, uint32_t(i), 0}); + } + if (lblTPC.isValid()) { + lookups.push_back({lblTPC, uint32_t(i), 1}); + } + if (mInterpolation.getTRDGIDsSuccess()[i].isIndexSet() && lblTPC.isValid()) { // track with TRD residuals: the true particle of the ITS-TPC part at the TRD + lookups.push_back({mc.label.isValid() ? mc.label : lblTPC, uint32_t(i), 2}); + } + if (lblTPC.isValid()) { // origin of the true particle of the ITS-TPC part + lookups.push_back({mc.label.isValid() ? mc.label : lblTPC, uint32_t(i), 3}); + } + } + // the reader loads the kinematics of a whole event (can be >100 MB): process event by event and release it right after + std::sort(lookups.begin(), lookups.end(), [](const Lookup& a, const Lookup& b) { + return a.lbl.getSourceID() != b.lbl.getSourceID() ? a.lbl.getSourceID() < b.lbl.getSourceID() : a.lbl.getEventID() < b.lbl.getEventID(); + }); + auto pdgToPID = [](int pdg) { + switch (std::abs(pdg)) { + case 11: + return o2::track::PID(o2::track::PID::Electron); + case 13: + return o2::track::PID(o2::track::PID::Muon); + case 321: + return o2::track::PID(o2::track::PID::Kaon); + case 2212: + return o2::track::PID(o2::track::PID::Proton); + case 1000010020: + return o2::track::PID(o2::track::PID::Deuteron); + case 1000010030: + return o2::track::PID(o2::track::PID::Triton); + case 1000020030: + return o2::track::PID(o2::track::PID::Helium3); + case 1000020040: + return o2::track::PID(o2::track::PID::Alpha); + default: + return o2::track::PID(o2::track::PID::Pion); + } + }; + auto refToPar = [&pdgToPID](const o2::TrackReference& ref, int charge, int pdg, bool sectorAlpha) { + std::array xyz{ref.X(), ref.Y(), ref.Z()}; + std::array pxyz{ref.Px(), ref.Py(), ref.Pz()}; + return o2::track::TrackPar(xyz, pxyz, charge, sectorAlpha, pdgToPID(pdg)); + }; + const auto matCorr = static_cast(mMatCorr); + auto prop = o2::base::Propagator::Instance(); + int curSrc = -1; + int curEv = -1; + for (const auto& lk : lookups) { + const auto& lbl = lk.lbl; + if (lbl.getSourceID() != curSrc || lbl.getEventID() != curEv) { + if (curSrc >= 0) { + mMCReader->releaseTracksForSourceAndEvent(curSrc, curEv); + } + curSrc = lbl.getSourceID(); + curEv = lbl.getEventID(); + } + const auto& trk = trkData[lk.idx]; + auto& mc = mTrackDataMC[lk.idx]; + const auto* mcTrk = mMCReader->getTrack(lbl); + int pdg = mcTrk ? mcTrk->GetPdgCode() : 0; + const auto* pPDG = mcTrk ? O2DatabasePDG::Instance()->GetParticle(pdg) : nullptr; + int charge = pPDG ? int(std::lround(pPDG->Charge() / 3.)) : 0; // TParticlePDG charge is in units of |e|/3 + auto refs = mMCReader->getTrackRefs(lbl.getSourceID(), lbl.getEventID(), lbl.getTrackID()); + if (lk.kind == 0) { // ITS outer: track reference of the ITS part nearest to TrackData::par + const o2::TrackReference* best = nullptr; + float bestD2 = 1e30f; + auto xyzReco = trk.par.getXYZGlo(); + for (const auto& ref : refs) { + if (ref.getDetectorId() != DetID::ITS) { + continue; + } + float dx = ref.X() - xyzReco.X(); + float dy = ref.Y() - xyzReco.Y(); + float dz = ref.Z() - xyzReco.Z(); + float d2 = dx * dx + dy * dy + dz * dz; + if (d2 < bestD2) { + bestD2 = d2; + best = &ref; + } + } + if (best && charge) { + auto par = refToPar(*best, charge, pdg, false); + if (par.rotateParam(trk.par.getAlpha()) && prop->PropagateToXBxByBz(par, trk.par.getX(), 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) { + mc.parITSOut = par; + mc.distITSRef = std::sqrt(bestD2); + mc.flags |= TrackDataMC::HasITSOut; + } + } + } else if (lk.kind == 1) { // TPC entrance: first TPC track reference in time of the TPC part + const o2::TrackReference* first = nullptr; + for (const auto& ref : refs) { + if (ref.getDetectorId() == DetID::TPC && (!first || ref.getTime() < first->getTime())) { + first = &ref; + } + } + if (first && charge) { + mc.parTPCIn = refToPar(*first, charge, pdg, true); + mc.flags |= TrackDataMC::HasTPCIn; + // association check (loopers, wrong leg, fake): truth at the innermost TPC cluster of the track vs that cluster + const auto& clRes = mInterpolation.getClusterResiduals(); + const UnbinnedResid* inner = nullptr; + for (int ic = trk.clIdx.getFirstEntry(); ic < trk.clIdx.getFirstEntry() + trk.clIdx.getEntries(); ++ic) { + if (clRes[ic].row < constants::MAXGLOBALPADROW && (!inner || clRes[ic].row < inner->row)) { + inner = &clRes[ic]; + } + } + if (inner) { + auto par = mc.parTPCIn; + float yCl = inner->y * param::MaxY / 0x7fff + inner->dy * param::MaxResid / 0x7fff; + float zCl = inner->z * param::MaxZ / 0x7fff + inner->dz * param::MaxResid / 0x7fff; + if (par.rotateParam(o2::math_utils::sector2Angle(inner->sec)) && prop->PropagateToXBxByBz(par, param::RowX[inner->row], 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) { + mc.distTPCRef = std::hypot(par.getY() - yCl, par.getZ() - zCl); + } + } + } + } else if (lk.kind == 3) { // origin: mother, production vertex and process + if (!mcTrk) { + continue; + } + mc.pdg = pdg; + if (mcTrk->isPrimary()) { + mc.flags |= TrackDataMC::IsPrimary; + } + mc.process = mcTrk->getProcess(); + mc.prodX = mcTrk->GetStartVertexCoordinatesX(); + mc.prodY = mcTrk->GetStartVertexCoordinatesY(); + mc.prodZ = mcTrk->GetStartVertexCoordinatesZ(); + mc.prodPx = mcTrk->GetStartVertexMomentumX(); + mc.prodPy = mcTrk->GetStartVertexMomentumY(); + mc.prodPz = mcTrk->GetStartVertexMomentumZ(); + int motherId = mcTrk->getMotherTrackId(); + if (motherId >= 0) { + mc.motherLabel = o2::MCCompLabel(motherId, lbl.getEventID(), lbl.getSourceID()); + const auto* mother = mMCReader->getTrack(lbl.getSourceID(), lbl.getEventID(), motherId); + mc.motherPdg = mother ? mother->GetPdgCode() : 0; + } + } else if (charge) { // TRD: entrance (first TRD track reference in time) and the true positions at the tracklet x of each layer + const o2::TrackReference* first = nullptr; + for (const auto& ref : refs) { + if (ref.getDetectorId() == DetID::TRD && (!first || ref.getTime() < first->getTime())) { + first = &ref; + } + } + if (!first) { + continue; + } + mc.parTRDIn = refToPar(*first, charge, pdg, true); + mc.flags |= TrackDataMC::HasTRDIn; + const auto& trkTRD = recoData.getITSTPCTRDTrack(mInterpolation.getTRDGIDsSuccess()[lk.idx]); // the TRD track of the stored TRD residuals + const auto trkltsCalib = recoData.getTRDCalibratedTracklets(); + const auto tracklets = recoData.getTRDTracklets(); + for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) { + int trkltIdx = trkTRD.getTrackletIndex(iLayer); + if (trkltIdx < 0) { + continue; + } + const auto& sp = trkltsCalib[trkltIdx]; // tracklet x, y, z in its sector frame, as used for the TRD residual + int sec = tracklets[trkltIdx].getDetector() / (o2::trd::constants::NLAYER * o2::trd::constants::NSTACK); + float alpha = o2::math_utils::sector2Angle(sec); + float cs = std::cos(alpha); + float sn = std::sin(alpha); + float gx = sp.getX() * cs - sp.getY() * sn; + float gy = sp.getX() * sn + sp.getY() * cs; + float gz = sp.getZ(); + const o2::TrackReference* best = nullptr; + float bestD2 = 1e30f; + for (const auto& ref : refs) { + if (ref.getDetectorId() != DetID::TRD) { + continue; + } + float dx = ref.X() - gx; + float dy = ref.Y() - gy; + float dz = ref.Z() - gz; + float d2 = dx * dx + dy * dy + dz * dz; + if (d2 < bestD2) { + bestD2 = d2; + best = &ref; + } + } + auto par = refToPar(*best, charge, pdg, false); + if (par.rotateParam(alpha) && prop->PropagateToXBxByBz(par, sp.getX(), 0.999f, o2::base::Propagator::MAX_STEP, matCorr)) { + mc.yTRD[iLayer] = par.getY(); + mc.zTRD[iLayer] = par.getZ(); + mc.trdLayerMask |= uint8_t(1) << iLayer; + } + } + } + } + if (curSrc >= 0) { + mMCReader->releaseTracksForSourceAndEvent(curSrc, curEv); + } + // stored tracks with the same mother (e.g. both daughters of a K0s): each points to the next one, cyclically + std::unordered_map> daughters; + for (uint32_t i = 0; i < mTrackDataMC.size(); ++i) { + if (mTrackDataMC[i].motherLabel.isSet()) { + daughters[mTrackDataMC[i].motherLabel.getTrackEventSourceID()].push_back(i); + } + } + for (const auto& [mother, idx] : daughters) { + for (size_t k = 0; idx.size() > 1 && k < idx.size(); ++k) { + mTrackDataMC[idx[k]].sisterIdx = idx[(k + 1) % idx.size()]; + } + } +} + void TPCInterpolationDPL::endOfStream(EndOfStreamContext& ec) { mInterpolation.finalize(); @@ -165,14 +428,13 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas dataRequest->setITSPerLayer(itsStag); std::vector outputs; - if (useMC) { - LOG(fatal) << "MC usage must be disabled for this workflow, since it is not yet implemented"; + dataRequest->requestTracks(srcVtx, false); + dataRequest->requestClusters(srcCls, false); + dataRequest->requestPrimaryVertices(false); + if (useMC) { // the MC truth needs only the labels of the ITS-TPC tracks and of their ITS and TPC parts + dataRequest->requestTracks(GTrackID::getSourcesMask("ITS,TPC,ITS-TPC"), true); } - dataRequest->requestTracks(srcVtx, useMC); - dataRequest->requestClusters(srcCls, useMC); - dataRequest->requestPrimaryVertices(useMC); - auto ggRequest = std::make_shared(false, // orbitResetTime true, // GRPECS=true true, // GRPLHCIF @@ -191,6 +453,9 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas if (debugOutput) { outputs.emplace_back("GLO", "TRKDATAEXT", 0, Lifetime::Timeframe); } + if (useMC && sendTrackData) { + outputs.emplace_back("GLO", "TRKDATAMC", 0, Lifetime::Timeframe); + } return DataProcessorSpec{ "tpc-track-interpolation", @@ -200,7 +465,8 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas Options{ {"matCorrType", VariantType::Int, 2, {"material correction type (definition in Propagator.h)"}}, {"sec-per-slot", VariantType::UInt32, 300u, {"number of seconds per calibration time slot (put 0 for infinite slot length)"}}, - {"process-seeds", VariantType::Bool, false, {"do not remove duplicates, e.g. for ITS-TPC-TRD track also process its seeding ITS-TPC part"}}}}; + {"process-seeds", VariantType::Bool, false, {"do not remove duplicates, e.g. for ITS-TPC-TRD track also process its seeding ITS-TPC part"}}, + {"mc-collision-context", VariantType::String, "collisioncontext.root", {"collision context used to access the MC kinematics and track references (MC only)"}}}}; } } // namespace tpc diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx index b3040d99bc4f2..a1c3bc7882c3d 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx @@ -189,6 +189,9 @@ void TPCResidualReader::run(ProcessingContext& pc) } for (int i = trkInfo.idxFirstResidual; i < trkInfo.idxFirstResidual + trkInfo.nResiduals; ++i) { const auto& residIn = mUnbinnedResiduals[i]; + if (residIn.isTgSlpClamped() || residIn.isPositionOnly()) { + continue; // scdcalib.clampTgSlp / keepClustersOnPropFail: tgSlp or dy, dz not usable for the binned voxel fit + } int sec = residIn.sec; auto& residVecOut = mResidualsSector[sec]; auto& statVecOut = mVoxStatsSector[sec]; diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx index e8b6aaa07eaba..d43c10ce9c74b 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx @@ -38,6 +38,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", VariantType::Bool, false, {"disable root-files input readers"}}, {"disable-root-output", VariantType::Bool, false, {"disable root-files output writers"}}, {"disable-mc", VariantType::Bool, false, {"disable MC propagation even if available"}}, + {"enable-mc", VariantType::Bool, false, {"store the MC truth of the track data (needs MC input and send-track-data)"}}, {"vtx-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources used for the vertex finding"}}, {"tracking-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use for track inter-/extrapolation"}}, {"tracking-sources-map-extraction", VariantType::String, std::string{GID::ALL}, {"can be subset of \"tracking-sources\""}}, @@ -103,8 +104,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); // write the configuration used for the workflow o2::conf::ConfigurableParam::writeINI("o2tpcinterpolation-workflow_configuration.ini"); - auto useMC = !configcontext.options().get("disable-mc"); - useMC = false; // force disabling MC as long as it is not implemented + // MC is opt-in: the residuals workflow is also run on data without passing disable-mc + auto useMC = configcontext.options().get("enable-mc") && !configcontext.options().get("disable-mc"); auto doStag = o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(configcontext); auto sendTrackData = configcontext.options().get("send-track-data"); auto debugOutput = configcontext.options().get("debug-output"); @@ -115,8 +116,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::tpc::getTPCResidualWriterSpec(sendTrackData, debugOutput)); } - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcClusters, srcVtx, srcVtx, useMC); - o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); // P-vertex is always needed + // MC labels only for the ITS-TPC tracks and their ITS and TPC parts (see getTPCInterpolationSpec) + GID::mask_t maskTracksMC = useMC ? GID::getSourcesMask("ITS,TPC,ITS-TPC") : GID::getSourcesMask(GID::NONE); + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcClusters, srcVtx, srcVtx, useMC, GID::getSourcesMask(GID::NONE), maskTracksMC); + o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, false); // P-vertex is always needed // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx index 20e37c3bcc3b4..d3593c1c50744 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-residual-aggregator.cxx @@ -32,6 +32,7 @@ void customize(std::vector& workflowOptions) {"output-type", VariantType::String, "unbinnedResid,trackParams", {"Comma separated list of outputs (without spaces). Valid strings: unbinnedResid, binnedResid, trackParams"}}, {"enable-track-input", VariantType::Bool, false, {"Whether to expect track data from interpolation workflow"}}, {"enable-ctp", VariantType::Bool, false, {"Subscribe to lumi info from CTP"}}, + {"enable-mc", VariantType::Bool, false, {"Whether to expect the MC truth of the track data from interpolation workflow (requires enable-track-input)"}}, {"disable-root-input", VariantType::Bool, false, {"disable root-files input readers"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); @@ -47,6 +48,14 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); auto trkInput = configcontext.options().get("enable-track-input"); auto ctpInput = configcontext.options().get("enable-ctp"); + auto mcInput = configcontext.options().get("enable-mc"); + if (mcInput && !trkInput) { + LOG(error) << "MC truth input requires the track input (enable-track-input), will be ignored"; + mcInput = false; + } + if (mcInput && !configcontext.options().get("disable-root-input")) { + LOG(fatal) << "MC truth input is only supported directly from the interpolation workflow (disable-root-input)"; + } bool writeUnbinnedResiduals = false; bool writeBinnedResiduals = false; @@ -78,7 +87,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) if (!configcontext.options().get("disable-root-input")) { specs.emplace_back(o2::tpc::getUnbinnedTPCResidualsReaderSpec(trkInput)); } - specs.emplace_back(getTPCResidualAggregatorSpec(trkInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData)); + specs.emplace_back(getTPCResidualAggregatorSpec(trkInput, ctpInput, writeUnbinnedResiduals, writeBinnedResiduals, writeTrackData, mcInput)); // CTP input if (ctpInput) { diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h index 00af697da3a9b..b2c5814a45f2e 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h @@ -45,11 +45,11 @@ struct ResidualsContainer { ResidualsContainer& operator=(const ResidualsContainer& src) = delete; ~ResidualsContainer(); - void init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, int autosave, int compression, long orbitResetTime); + void init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, bool wTrackDataMC, int autosave, int compression, long orbitResetTime); void fillStatisticsBranches(); uint64_t getNEntries() const { return nResidualsTotal; } - void fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const o2::ctp::LumiInfo* lumiInput); + void fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const gsl::span* trkDataMCIn, const o2::ctp::LumiInfo* lumiInput); void merge(ResidualsContainer* prev); void print(); void writeToFile(bool closeFileAfterwards); @@ -66,6 +66,7 @@ struct ResidualsContainer { std::vector unbinnedRes, *unbinnedResPtr{&unbinnedRes}; ///< unbinned residuals which are sent to the aggregator std::vector detInfoUnbRes, *detInfoUnbResPtr{&detInfoUnbRes}; ///< detector info associated to unbinned residuals which are sent to the aggregator std::vector trkData, *trkDataPtr{&trkData}; ///< track data and cluster ranges + std::vector trkDataMC, *trkDataMCPtr{&trkDataMC}; ///< MC truth for the track data (MC only) std::vector trackInfo, *trackInfoPtr{&trackInfo}; ///< allows to obtain track type for each unbinned residual downstream o2::ctp::LumiInfo lumiTF; ///< for each processed TF we store the lumi information in the tree of unbinned residuals uint64_t timeMS; ///< for each processed TF we store its absolute time in ms in the tree of unbinned residuals @@ -83,6 +84,7 @@ struct ResidualsContainer { bool writeBinnedResid{false}; ///< flag, whether binned residuals should be written out bool writeUnbinnedResiduals{false}; ///< flag, whether unbinned residuals should be written out bool writeTrackData{false}; ///< flag, whether full seeding track information should be written out + bool writeTrackDataMC{false}; ///< flag, whether the MC truth of the seeding tracks should be written out int autosaveInterval{0}; ///< if > 0, then the output written to file for every n-th TF // additional info @@ -94,7 +96,7 @@ struct ResidualsContainer { float TPCVDriftRef{-1.}; ///< TPC nominal drift speed in cm/microseconds float TPCDriftTimeOffsetRef{0.}; ///< TPC nominal (e.g. at the start of run) drift time bias in cm/mus - ClassDefNV(ResidualsContainer, 5); + ClassDefNV(ResidualsContainer, 6); }; class ResidualAggregator final : public o2::calibration::TimeSlotCalibration @@ -122,6 +124,7 @@ class ResidualAggregator final : public o2::calibration::TimeSlotCalibration0 then the output is written to a file for every n-th TF int mCompressionSetting{101}; ///< single integer defining the ROOT compression algorithm and level (see TFile doc for details) size_t mMinEntries; ///< the minimum number of residuals required for the map creation (per voxel) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h index 6ef3839991d04..7bba240552215 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h @@ -49,6 +49,7 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelper= param::MaxTgSlp saturated (tgSlp = +-0x7fff, see UnbinnedResid::isTgSlpClamped) instead of dropping them: the cut is on the reference track's direction, so dropping selects on the reference's error + bool keepClustersOnPropFail{false}; ///< if the reference track propagation fails (maxSnp, rotation), keep the track: TPC clusters without a reference are stored position-only (y, z = cluster; dy = dz = 0; tgSlp = UnbinnedResid::TgSlpPositionOnly) instead of dropping the whole track, which selects on the reference's error float maxStep{2.f}; ///< maximum step for propagation bool debugTRDTOF{false}; ///< if true, ITS-TPC-TRD-TOF tracks and their seeding ITS-TPC-TRD track will both be interpolated and their residuals stored diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index 91726aa9941fa..38b7fd3784462 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -24,6 +24,7 @@ #include "ReconstructionDataFormats/TrackTPCITS.h" #include "ReconstructionDataFormats/MatchInfoTOF.h" #include "ReconstructionDataFormats/GlobalTrackID.h" +#include "SimulationDataFormat/MCCompLabel.h" #include "DataFormatsITSMFT/Cluster.h" #include "DataFormatsITSMFT/TrkClusRef.h" #include "DataFormatsITSMFT/TopologyDictionary.h" @@ -97,6 +98,10 @@ struct UnbinnedResid { /// true if tgSlp was saturated at +-param::MaxTgSlp (scdcalib.clampTgSlp): unclamped values have |tgSlp| <= 0x7fff - 1 bool isTgSlpClamped() const { return tgSlp == 0x7fff || tgSlp == -0x7fff; } + /// tgSlp marker of a position-only TPC cluster (scdcalib.keepClustersOnPropFail): no reference track at this cluster, + /// y and z are the cluster position, dy = dz = 0. Not reachable by the tgSlp packing (|tgSlp| <= 0x7fff) + static constexpr short TgSlpPositionOnly = -0x8000; + bool isPositionOnly() const { return tgSlp == TgSlpPositionOnly; } bool isTPC() const { return row < constants::MAXGLOBALPADROW; } bool isTRD() const { return row >= 160 && row < 166; } bool isTOF() const { return row == 170; } @@ -240,6 +245,44 @@ struct TrackData { ClassDefNV(TrackData, 12); }; +/// MC truth for a TrackData entry (stored only for MC, aligned 1:1 with the TrackData vector) +struct TrackDataMC { + enum Flags : uint8_t { HasITSOut = 0x1, ///< parITSOut is filled + HasTPCIn = 0x2, ///< parTPCIn is filled + FakeITSTPC = 0x4, ///< ITS and TPC parts of the track have different MC labels + HasTRDIn = 0x8, ///< parTRDIn is filled + IsPrimary = 0x10 }; ///< the particle of the ITS-TPC part is a primary (MCTrack::isPrimary) + o2::MCCompLabel label{}; ///< MC label of the ITS-TPC part of the seeding track + o2::MCCompLabel labelITS{}; ///< MC label of its ITS part + o2::MCCompLabel labelTPC{}; ///< MC label of its TPC part + o2::track::TrackPar parITSOut{}; ///< truth at x and alpha of TrackData::par, from the nearest ITS track reference (propagated with the material correction, true mass) + o2::track::TrackPar parTPCIn{}; ///< truth at the first TPC track reference (sector frame) + float distITSRef{-1.f}; ///< 3D distance between the ITS track reference used and TrackData::par in cm + float distTPCRef{-1.f}; ///< distance (y,z) between parTPCIn propagated to the innermost TPC cluster of the track and that cluster in cm (large: wrong leg, looper, fake) + o2::track::TrackPar parTRDIn{}; ///< truth at the first TRD track reference (sector frame), TRD-matched seeds only + float yTRD[6] = {}; ///< truth y at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask + float zTRD[6] = {}; ///< truth z at the x of the TRD tracklet of each layer (tracklet sector frame), see trdLayerMask + uint8_t trdLayerMask{0}; ///< bit i set: yTRD[i], zTRD[i] filled + int pdg{0}; ///< PDG code of the particle of the ITS-TPC part (as for the origin and TRD fields) + o2::MCCompLabel motherLabel{}; ///< MC label of the mother of the particle of the ITS-TPC part (for primaries the generator-level parent) + int motherPdg{0}; ///< PDG code of that mother (0: none) + float prodX{0.f}; ///< production vertex x of the particle of the ITS-TPC part (global, cm) + float prodY{0.f}; ///< production vertex y (global, cm) + float prodZ{0.f}; ///< production vertex z (global, cm) + float prodPx{0.f}; ///< momentum at production x (global, GeV/c), e.g. to compare a track propagated to the vertex + float prodPy{0.f}; ///< momentum at production y (global, GeV/c) + float prodPz{0.f}; ///< momentum at production z (global, GeV/c) + int sisterIdx{-1}; ///< index in the TrackData vector of another stored track with the same mother (cycling through all of them if more than two), -1: none + uint8_t process{0}; ///< production process of the particle of the ITS-TPC part (TMCProcess) + uint8_t flags{0}; + bool hasITSOut() const { return flags & HasITSOut; } + bool hasTPCIn() const { return flags & HasTPCIn; } + bool isFakeITSTPC() const { return flags & FakeITSTPC; } + bool hasTRDIn() const { return flags & HasTRDIn; } + bool isPrimary() const { return flags & IsPrimary; } + ClassDefNV(TrackDataMC, 3); +}; + /// \class TrackInterpolation /// This class is retrieving the TPC space point residuals by interpolating ITS/TRD/TOF tracks. /// The residuals are stored in the specified vectors of TPCClusterResiduals @@ -420,6 +463,8 @@ class TrackInterpolation std::vector& getTrackDataCompact() { return mTrackDataCompact; } std::vector& getTrackDataExtended() { return mTrackDataExtended; } std::vector& getReferenceTracks() { return mTrackData; } + /// ITS-TPC-TRD track whose tracklets gave the TRD residuals of each stored track (not set if none), aligned with getReferenceTracks() + const std::vector& getTRDGIDsSuccess() const { return mTRDGIDsSuccess; } void setLane(int lID, int nL) { @@ -487,6 +532,7 @@ class TrackInterpolation // cache std::array mCache{{}}; ///< caching positions, covariances and angles for track extrapolations and interpolation std::vector mGIDsSuccess; ///< keep track of the GIDs which could be processed successfully + std::vector mTRDGIDsSuccess; ///< ITS-TPC-TRD track used for the TRD residuals of each stored track (not set if none) TrackValidationData mTrackValidation; @@ -500,6 +546,8 @@ class TrackInterpolation size_t mNRejRefit = 0; size_t mNRejProp = 0; size_t mNRejLoop = 0; + size_t mNPosOnlyTracks = 0; ///< tracks kept with position-only clusters after a propagation failure (keepClustersOnPropFail) + size_t mNPosOnlyClusters = 0; ///< position-only TPC clusters stored (keepClustersOnPropFail) ClassDefNV(TrackInterpolation, 1); }; diff --git a/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C b/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C index 3bf01f21f1dce..ee5c54f32acf5 100644 --- a/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C +++ b/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C @@ -320,9 +320,9 @@ void staticMapCreator(std::string fileInput = "files.txt", if (useResidualsForVd && residualsVd.size() < 10'000'000UL) { residualsVd.push_back(residIn); } - if (residIn.isTgSlpClamped()) { + if (residIn.isTgSlpClamped() || residIn.isPositionOnly()) { // scdcalib.clampTgSlp: tgSlp saturated -- the voxel fit (dX from dY vs tan(phi)) and the map correction below - // use it, so keep this residual out of the binned residuals + // use it, so keep this residual out of the binned residuals; scdcalib.keepClustersOnPropFail: no reference (dy = dz = 0) continue; } int sec = residIn.sec; diff --git a/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx b/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx index c5594ddc40b02..e8423957990c7 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/ResidualAggregator.cxx @@ -77,6 +77,7 @@ ResidualsContainer::ResidualsContainer(ResidualsContainer&& rhs) unbinnedRes = std::move(rhs.unbinnedRes); trackInfo = std::move(rhs.trackInfo); trkData = std::move(rhs.trkData); + trkDataMC = std::move(rhs.trkDataMC); orbitReset = rhs.orbitReset; firstTForbit = rhs.firstTForbit; firstSeenTF = rhs.firstSeenTF; @@ -84,13 +85,14 @@ ResidualsContainer::ResidualsContainer(ResidualsContainer&& rhs) nResidualsTotal = rhs.nResidualsTotal; } -void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, int autosave, int compression, long orbitResetTime) +void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string outputDir, bool wFile, bool wBinnedResid, bool wUnbinnedResid, bool wTrackData, bool wTrackDataMC, int autosave, int compression, long orbitResetTime) { trackResiduals = residualsEngine; writeToRootFile = wFile; writeBinnedResid = wBinnedResid; writeUnbinnedResiduals = wUnbinnedResid; writeTrackData = wTrackData; + writeTrackDataMC = wTrackData && wTrackDataMC; autosaveInterval = autosave; orbitReset = orbitResetTime; if (writeToRootFile) { @@ -129,6 +131,9 @@ void ResidualsContainer::init(const TrackResiduals* residualsEngine, std::string if (writeTrackData) { treeOutTrackData = std::make_unique("trackData", "Track information incl cluster range ref"); treeOutTrackData->Branch("trk", &trkDataPtr); + if (writeTrackDataMC) { + treeOutTrackData->Branch("trkMC", &trkDataMCPtr); + } } if (writeBinnedResid) { treeOutResiduals = std::make_unique("resid", "TPC binned residuals"); @@ -171,7 +176,7 @@ void ResidualsContainer::fillStatisticsBranches() } } -void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const o2::ctp::LumiInfo* lumiInput) +void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::span resid, const gsl::span detInfoRes, const gsl::span trkRefsIn, const gsl::span* trkDataIn, const gsl::span* trkDataMCIn, const o2::ctp::LumiInfo* lumiInput) { // receives large vector of unbinned residuals and fills the sector-wise vectors // with binned residuals and statistics @@ -197,9 +202,9 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp if (!writeBinnedResid) { continue; } - if (residIn.isTgSlpClamped()) { + if (residIn.isTgSlpClamped() || residIn.isPositionOnly()) { // scdcalib.clampTgSlp: kept in the unbinned output, but its tgSlp is saturated and the voxel fit uses tgSlp (dX from - // dY vs tan(phi)), so it must not enter the binned residuals + // dY vs tan(phi)), so it must not enter the binned residuals; scdcalib.keepClustersOnPropFail: no reference, dy = dz = 0 continue; } int sec = residIn.sec; @@ -244,8 +249,12 @@ void ResidualsContainer::fill(const o2::dataformats::TFIDInfo& ti, const gsl::sp for (const auto& trkIn : *trkDataIn) { trkData.push_back(trkIn); } + if (writeTrackDataMC && trkDataMCIn) { + trkDataMC.assign(trkDataMCIn->begin(), trkDataMCIn->end()); + } treeOutTrackData->Fill(); trkData.clear(); + trkDataMC.clear(); } if (writeUnbinnedResiduals) { if (lumiInput) { @@ -338,6 +347,9 @@ void ResidualsContainer::merge(ResidualsContainer* prev) if (writeTrackData) { prev->treeOutTrackData->SetBranchAddress("trk", &trkDataPtr); + if (writeTrackDataMC) { + prev->treeOutTrackData->SetBranchAddress("trkMC", &trkDataMCPtr); + } for (int i = 0; i < treeOutTrackData->GetEntries(); ++i) { treeOutTrackData->GetEntry(i); prev->treeOutTrackData->Fill(); @@ -456,7 +468,7 @@ Slot& ResidualAggregator::emplaceNewSlot(bool front, TFType tStart, TFType tEnd) auto& cont = getSlots(); auto& slot = front ? cont.emplace_front(tStart, tEnd) : cont.emplace_back(tStart, tEnd); slot.setContainer(std::make_unique()); - slot.getContainer()->init(&mTrackResiduals, mOutputDir, mWriteOutput, mWriteBinnedResiduals, mWriteUnbinnedResiduals, mWriteTrackData, mAutosaveInterval, mCompressionSetting, mOrbitResetTime); + slot.getContainer()->init(&mTrackResiduals, mOutputDir, mWriteOutput, mWriteBinnedResiduals, mWriteUnbinnedResiduals, mWriteTrackData, mWriteTrackDataMC, mAutosaveInterval, mCompressionSetting, mOrbitResetTime); std::chrono::duration emplaceDuration = std::chrono::high_resolution_clock::now() - emplaceStartTime; LOGP(info, "Emplacing new calibration slot took: {} ms", std::chrono::duration_cast(emplaceDuration).count()); return slot; diff --git a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h index e77610acb8e7e..b6db3abfecc90 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h +++ b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h @@ -21,6 +21,8 @@ #pragma link C++ class std::vector < o2::tpc::TrackDataCompact> + ; #pragma link C++ class o2::tpc::TrackData + ; #pragma link C++ class std::vector < o2::tpc::TrackData> + ; +#pragma link C++ class o2::tpc::TrackDataMC + ; +#pragma link C++ class std::vector < o2::tpc::TrackDataMC> + ; #pragma link C++ class o2::tpc::TrackDataExtended + ; #pragma link C++ class std::vector < o2::tpc::TrackDataExtended> + ; #pragma link C++ class o2::tpc::TPCClusterResiduals + ; diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx index 86fca56e4f8eb..17237473e85ae 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx @@ -465,12 +465,26 @@ void TrackInterpolation::process() } LOGP(info, "Could process {} tracks successfully ({} rejected in refits, {} in propagation, {} as loopers), {} residuals were rejected, {} accepted", mTrackData.size(), mNRejRefit, mNRejProp, mNRejLoop, mRejectedResiduals, mClRes.size()); + if (mParams->keepClustersOnPropFail) { + LOGP(info, "keepClustersOnPropFail: {} tracks kept after a propagation failure, {} position-only TPC clusters stored", mNPosOnlyTracks, mNPosOnlyClusters); + } + mNPosOnlyTracks = 0; + mNPosOnlyClusters = 0; mRejectedResiduals = 0; mNRejRefit = 0; mNRejProp = 0; mNRejLoop = 0; } +namespace +{ +/// TPC cluster stored without a reference track (scdcalib.keepClustersOnPropFail), see UnbinnedResid::isPositionOnly +struct PositionOnlyCluster { + float y, z; + unsigned char sec, row, flags; +}; +} // namespace + void TrackInterpolation::interpolateTrack(int iSeed) { LOGP(debug, "Starting track interpolation for GID {}", mGIDs[iSeed].asString()); @@ -514,6 +528,11 @@ void TrackInterpolation::interpolateTrack(int iSeed) // store the TPC cluster positions in the cache, as well as dedx info std::array, constants::MAXGLOBALPADROW> mCacheDEDX{}; std::array multBins{}; + // keepClustersOnPropFail: a row gets a residual only if both the outward (ITS) and the inward (TRD/TOF) propagation reached + // it; the other TPC clusters are stored position-only. allLost: the outward pass or the outer anchor failed. + bool allLost = false; + std::array refOut{}, refIn{}; + std::vector posOnly; for (int iCl = trkTPC.getNClusterReferences(); iCl--;) { uint8_t sector, row; uint32_t clusterIndexInRow; @@ -545,16 +564,16 @@ void TrackInterpolation::interpolateTrack(int iSeed) if (!mCache[iRow].clAvailable) { continue; } - if (!trkWork.rotate(mCache[iRow].clAngle)) { - LOG(debug) << "Failed to rotate track during first extrapolation"; - mNRejProp++; - return; - } - if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { + if (!trkWork.rotate(mCache[iRow].clAngle) || !propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed on first extrapolation"; - mNRejProp++; - return; + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + allLost = true; // no outer anchor can be reached: every TPC cluster is stored position-only + break; } + refOut[iRow] = true; mCache[iRow].y[ExtOut] = trkWork.getY(); mCache[iRow].z[ExtOut] = trkWork.getZ(); mCache[iRow].sy2[ExtOut] = trkWork.getSigmaY2(); @@ -565,7 +584,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } // start from outermost cluster with outer refit and back propagation - if (gidTable[GTrackID::TOF].isIndexSet()) { + if (!allLost && gidTable[GTrackID::TOF].isIndexSet()) { LOG(debug) << "TOF point available"; const auto& clTOF = mRecoCont->getTOFClusters()[gidTable[GTrackID::TOF]]; if (mDumpTrackPoints) { @@ -574,31 +593,23 @@ void TrackInterpolation::interpolateTrack(int iSeed) } const int clTOFSec = clTOF.getCount(); const float clTOFAlpha = o2::math_utils::sector2Angle(clTOFSec); - if (!trkWork.rotate(clTOFAlpha)) { - LOG(debug) << "Failed to rotate into TOF cluster sector frame"; - mNRejProp++; - return; - } float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()}; if (!clTOF.isInNominalSector()) { o2::tof::Geo::alignedToNominalSector(clTOFxyz, clTOFSec); // go from the aligned to nominal sector frame } std::array clTOFYZ{clTOFxyz[1], clTOFxyz[2]}; std::array clTOFCov{mParams->sigYZ2TOF, 0.f, mParams->sigYZ2TOF}; // assume no correlation between y and z and equal cluster error sigma^2 = (3cm)^2 / 12 - if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) { - LOG(debug) << "Failed final propagation to TOF radius"; - mNRejProp++; - return; - } // TODO: check if reset of covariance matrix is needed here (or, in case TOF point is not available at outermost TRD layer) - if (!trkWork.update(clTOFYZ, clTOFCov)) { - LOG(debug) << "Failed to update extrapolated ITS track with TOF cluster"; - // LOGF(info, "trkWork.y=%f, cl.y=%f, trkWork.z=%f, cl.z=%f", trkWork.getY(), clTOFYZ[0], trkWork.getZ(), clTOFYZ[1]); - mNRejProp++; - return; + if (!trkWork.rotate(clTOFAlpha) || !propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr) || !trkWork.update(clTOFYZ, clTOFCov)) { + LOG(debug) << "Failed to rotate/propagate/update the extrapolated ITS track at the TOF cluster"; + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + allLost = true; } } - if (gidTable[GTrackID::TRD].isIndexSet()) { + if (!allLost && gidTable[GTrackID::TRD].isIndexSet()) { LOG(debug) << "TRD available"; const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTable[GTrackID::ITSTPCTRD]); if (mDumpTrackPoints) { @@ -611,13 +622,16 @@ void TrackInterpolation::interpolateTrack(int iSeed) if (res == -1) { // no TRD tracklet in this layer continue; } - if (res < -1) { // failed to reach this layer - return; - } - if (!trkWork.update(trkltTRDYZ, trkltTRDCov)) { - LOG(debug) << "Failed to update track at TRD layer " << iLayer; - mNRejProp++; - return; + if (res < -1 || !trkWork.update(trkltTRDYZ, trkltTRDCov)) { // failed to reach this layer or to update + LOG(debug) << "Failed to reach or update the track at TRD layer " << iLayer; + if (!mParams->keepClustersOnPropFail) { + if (res >= -1) { + mNRejProp++; // unchanged: only the update failure was counted + } + return; + } + allLost = true; + break; } } } @@ -629,7 +643,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) // go back through the TPC and store updated track positions bool outerParamStored = false; - for (int iRow = param::NPadRows; iRow--;) { + for (int iRow = param::NPadRows; !allLost && iRow--;) { if (!mCache[iRow].clAvailable) { continue; } @@ -642,17 +656,15 @@ void TrackInterpolation::interpolateTrack(int iSeed) trackData.par = trkWork; outerParamStored = true; } - if (!trkWork.rotate(mCache[iRow].clAngle)) { - LOG(debug) << "Failed to rotate track during back propagation"; - mNRejProp++; - return; - } - if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { + if (!trkWork.rotate(mCache[iRow].clAngle) || !propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed on back propagation"; - // printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha); - mNRejProp++; - return; + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + break; // this row and all inner ones have no inward reference: stored position-only } + refIn[iRow] = true; mCache[iRow].y[ExtIn] = trkWork.getY(); mCache[iRow].z[ExtIn] = trkWork.getZ(); mCache[iRow].sy2[ExtIn] = trkWork.getSigmaY2(); @@ -668,6 +680,11 @@ void TrackInterpolation::interpolateTrack(int iSeed) ++deltaRow; continue; } + if (!refOut[iRow] || !refIn[iRow]) { // keepClustersOnPropFail only: no reference at this row + posOnly.push_back({mCache[iRow].clY, mCache[iRow].clZ, mCache[iRow].clSec, (unsigned char)iRow, mCache[iRow].clFlags}); + ++deltaRow; + continue; + } float wTotY = 1.f / mCache[iRow].sy2[ExtOut] + 1.f / mCache[iRow].sy2[ExtIn]; float wTotZ = 1.f / mCache[iRow].sz2[ExtOut] + 1.f / mCache[iRow].sz2[ExtIn]; mCache[iRow].y[Int] = (mCache[iRow].y[ExtOut] / mCache[iRow].sy2[ExtOut] + mCache[iRow].y[ExtIn] / mCache[iRow].sy2[ExtIn]) / wTotY; @@ -713,12 +730,30 @@ void TrackInterpolation::interpolateTrack(int iSeed) mTrackValidation.clear(); // for refitted track parameters and flagging rejected clusters bool stored = false; - trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals, true); + // keepClustersOnPropFail: a track without any interpolated residual has nothing to validate + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : ((mParams->keepClustersOnPropFail && clusterResiduals.empty()) ? int8_t(0x1) : validateTrack(trackData, mTrackValidation, clusterResiduals, true)); if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0; int iRow = 0; + // keepClustersOnPropFail: store the position-only clusters in row order between the residuals + size_t iPosOnly = 0; + auto flushPosOnly = [&](int rowLimit) { + for (; iPosOnly < posOnly.size() && posOnly[iPosOnly].row < rowLimit; ++iPosOnly) { + const auto& pc = posOnly[iPosOnly]; + if (std::abs(pc.y) < param::MaxY && std::abs(pc.z) < param::MaxZ) { + mClRes.emplace_back(0.f, 0.f, 0.f, pc.y, pc.z, pc.row, pc.sec, pc.flags, false); + mClRes.back().tgSlp = UnbinnedResid::TgSlpPositionOnly; + mDetInfoRes.emplace_back().setTPC(mCacheDEDX[pc.row].first, mCacheDEDX[pc.row].second); // qtot, qmax + ++nClValidated; + ++mNPosOnlyClusters; + } else { + ++mRejectedResiduals; + } + } + }; for (unsigned int iCl = 0; iCl < clusterResiduals.size(); ++iCl) { iRow += clusterResiduals[iCl].dRow; + flushPosOnly(iRow); const auto rej = trackData.filterFlag < 0 ? false : mTrackValidation.points[iCl].flagRej; if (rej && !mParams->keepRejectedResiduals) { // skip masked cluster residual continue; @@ -741,6 +776,10 @@ void TrackInterpolation::interpolateTrack(int iSeed) ++mRejectedResiduals; } } + flushPosOnly(constants::MAXGLOBALPADROW); + if (!posOnly.empty()) { + ++mNPosOnlyTracks; + } trackData.clIdx.setEntries(nClValidated); // store multiplicity info @@ -768,10 +807,12 @@ void TrackInterpolation::interpolateTrack(int iSeed) } bool stopPropagation = !mExtDetResid; + GTrackID gidTRDUsed{}; if (!stopPropagation) { // do we have TRD residuals to add? trkWork = trkOuter; - if (gidTable[GTrackID::TRD].isIndexSet()) { + if (!allLost && gidTable[GTrackID::TRD].isIndexSet()) { // allLost: trkOuter is not a valid outer param + gidTRDUsed = gidTable[GTrackID::ITSTPCTRD]; const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTable[GTrackID::ITSTPCTRD]); for (int iLayer = 0; iLayer < o2::trd::constants::NLAYER; iLayer++) { std::array trkltTRDYZ{}; @@ -796,7 +837,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } // do we have TOF residual to add? - while (gidTable[GTrackID::TOF].isIndexSet() && !stopPropagation) { + while (!allLost && gidTable[GTrackID::TOF].isIndexSet() && !stopPropagation) { const auto& clTOF = mRecoCont->getTOFClusters()[gidTable[GTrackID::TOF]]; float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()}; if (!clTOF.isInNominalSector()) { @@ -887,6 +928,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } mGIDsSuccess.push_back(mGIDs[iSeed]); + mTRDGIDsSuccess.push_back(gidTRDUsed); mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag); mTrackData.push_back(std::move(trackData)); stored = true; @@ -996,6 +1038,8 @@ void TrackInterpolation::extrapolateTrack(int iSeed) uint8_t clRowPrev = constants::MAXGLOBALPADROW; // used to identify and skip split clusters on the same pad row std::array, constants::MAXGLOBALPADROW> mCacheDEDX{}; std::array multBins{}; + bool refLost = false; // keepClustersOnPropFail: the ITS extrapolation failed at an earlier cluster + std::vector posOnly; // keepClustersOnPropFail: clusters after the failure, rows ascending for (int iCl = trkTPC.getNClusterReferences(); iCl--;) { uint8_t sector, row; uint32_t clusterIndexInRow; @@ -1017,29 +1061,31 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } float x = 0, y = 0, z = 0; mFastTransform->TransformIdeal(sector, row, cl.getPad(), cl.getTime(), x, y, z, clusterTimeBinOffset); - if (!trkWork.rotate(o2::math_utils::sector2Angle(sector))) { - mNRejProp++; - return; - } - if (!propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr)) { - mNRejProp++; - return; - } - - const auto dY = y - trkWork.getY(); - const auto dZ = z - trkWork.getZ(); - const auto ty = trkWork.getY(); - const auto tz = trkWork.getZ(); - const auto snp = trkWork.getSnp(); - const auto sec = sector; unsigned char flags = cl.getFlags(); if (mTPCShClassMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) { flags |= o2::gpu::GPUTPCGMMergedTrackHit::flagShared; } - clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev, flags); + if (!refLost && !(trkWork.rotate(o2::math_utils::sector2Angle(sector)) && propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr))) { + if (!mParams->keepClustersOnPropFail) { + mNRejProp++; + return; + } + refLost = true; // scdcalib.keepClustersOnPropFail: this and all further clusters are stored position-only + } mCacheDEDX[row].first = cl.getQtot(); mCacheDEDX[row].second = cl.getQmax(); - rowPrev = row; + if (refLost) { + posOnly.push_back({y, z, sector, row, flags}); + } else { + const auto dY = y - trkWork.getY(); + const auto dZ = z - trkWork.getZ(); + const auto ty = trkWork.getY(); + const auto tz = trkWork.getZ(); + const auto snp = trkWork.getSnp(); + const auto sec = sector; + clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev, flags); + rowPrev = row; + } int imb = int(cl.getTime() * mNTPCOccBinLengthInv); if (imb < mTPCParam->occupancyMapSize) { multBins[row] = 1 + std::max(0, imb); @@ -1065,12 +1111,30 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } bool stored = false; - trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals, false); + // keepClustersOnPropFail: a track that lost its reference before the first cluster has no residual to validate + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : ((mParams->keepClustersOnPropFail && clusterResiduals.empty()) ? int8_t(0x1) : validateTrack(trackData, mTrackValidation, clusterResiduals, false)); if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0, iRow = 0; unsigned int iCl = 0; + // keepClustersOnPropFail: store the position-only clusters in row order between the residuals + size_t iPosOnly = 0; + auto flushPosOnly = [&](int rowLimit) { + for (; iPosOnly < posOnly.size() && posOnly[iPosOnly].row < rowLimit; ++iPosOnly) { + const auto& pc = posOnly[iPosOnly]; + if (std::abs(pc.y) < param::MaxY && std::abs(pc.z) < param::MaxZ) { + mClRes.emplace_back(0.f, 0.f, 0.f, pc.y, pc.z, pc.row, pc.sec, pc.flags, false); + mClRes.back().tgSlp = UnbinnedResid::TgSlpPositionOnly; + mDetInfoRes.emplace_back().setTPC(mCacheDEDX[pc.row].first, mCacheDEDX[pc.row].second); // qtot, qmax + ++nClValidated; + ++mNPosOnlyClusters; + } else { + ++mRejectedResiduals; + } + } + }; for (iCl = 0; iCl < clusterResiduals.size(); ++iCl) { iRow += clusterResiduals[iCl].dRow; + flushPosOnly(iRow); if (iRow >= param::NPadRows) { // RS why do we need this? continue; } @@ -1095,6 +1159,10 @@ void TrackInterpolation::extrapolateTrack(int iSeed) ++mRejectedResiduals; } } + flushPosOnly(constants::MAXGLOBALPADROW); + if (!posOnly.empty()) { + ++mNPosOnlyTracks; + } trackData.clIdx.setEntries(nClValidated); // store multiplicity info @@ -1122,12 +1190,14 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } bool stopPropagation = !mExtDetResid; + GTrackID gidTRDUsed{}; if (!stopPropagation) { // do we have TRD residuals to add? int iSeedFull = mParentID[iSeed] == -1 ? iSeed : mParentID[iSeed]; auto gidFull = mGIDs[iSeedFull]; const auto& gidTableFull = mGIDtables[iSeedFull]; - if (gidTableFull[GTrackID::TRD].isIndexSet()) { + if (!refLost && gidTableFull[GTrackID::TRD].isIndexSet()) { // refLost: trkWork did not reach the TPC outer end + gidTRDUsed = gidTableFull[GTrackID::ITSTPCTRD]; const auto& trkTRD = mRecoCont->getITSTPCTRDTrack(gidTableFull[GTrackID::ITSTPCTRD]); trackData.nTrkltsTRD = trkTRD.getNtracklets(); trackData.chi2TRD = trkTRD.getChi2(); @@ -1156,7 +1226,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) // do we have TOF residual to add? trackData.clAvailTOF = 0; - while (gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) { + while (!refLost && gidTableFull[GTrackID::TOF].isIndexSet() && !stopPropagation) { const auto& tofMatch = mRecoCont->getTOFMatch(gidFull); ULong64_t bclongtof = (tofMatch.getSignal() - 10000) * o2::tof::Geo::BC_TIME_INPS_INV; double t0forTOF = tofMatch.getFT0Best(); // setting t0 for TOF @@ -1254,6 +1324,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) mTrackData.push_back(std::move(trackData)); stored = true; mGIDsSuccess.push_back(mGIDs[iSeed]); + mTRDGIDsSuccess.push_back(gidTRDUsed); mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag); if (mDumpTrackPoints) { (*trackDataExtended).clIdx.setEntries(nClValidated); @@ -1629,6 +1700,7 @@ void TrackInterpolation::reset() mClRes.clear(); mDetInfoRes.clear(); mGIDsSuccess.clear(); + mTRDGIDsSuccess.clear(); for (auto& vec : mTrackIndices) { vec.clear(); } diff --git a/Steer/include/Steer/MCKinematicsReader.h b/Steer/include/Steer/MCKinematicsReader.h index ae5ccf6615c56..793711c61de87 100644 --- a/Steer/include/Steer/MCKinematicsReader.h +++ b/Steer/include/Steer/MCKinematicsReader.h @@ -87,7 +87,7 @@ class MCKinematicsReader /// variant returning all tracks for source and event at once std::vector const& getTracks(int source, int event) const; - /// API to ask releasing tracks (freeing memory) for source + event + /// API to ask releasing tracks and track references (freeing memory) for source + event void releaseTracksForSourceAndEvent(int source, int event); /// variant returning all tracks for an event id (source = 0) at once @@ -129,7 +129,8 @@ class MCKinematicsReader void initTracksForSource(int source) const; void loadTracksForSourceAndEvent(int source, int eventID) const; void loadHeadersForSource(int source) const; - void loadTrackRefsForSource(int source) const; + void initTrackRefsForSource(int source) const; + void loadTrackRefsForSourceAndEvent(int source, int event) const; void initIndexedTrackRefs(std::vector& refs, o2::dataformats::MCTruthContainer& indexedrefs) const; DigitizationContext const* mDigitizationContext = nullptr; @@ -142,6 +143,7 @@ class MCKinematicsReader mutable std::vector*>> mTracks; // the in-memory track container mutable std::vector> mHeaders; // the in-memory header container mutable std::vector>> mIndexedTrackRefs; // the in-memory track ref container + mutable std::vector> mTrackRefsLoaded; // whether the track refs of a source/event are in memory bool mInitialized = false; // whether initialized }; @@ -206,11 +208,14 @@ inline gsl::span MCKinematicsReader::getTrackRefs(int source } auto& perEvent = mIndexedTrackRefs[source]; if (perEvent.size() == 0) { - loadTrackRefsForSource(source); + initTrackRefsForSource(source); } if (static_cast(event) >= perEvent.size()) { return {}; } + if (!mTrackRefsLoaded[source][event]) { + loadTrackRefsForSourceAndEvent(source, event); + } return perEvent[event].getLabels(track); } @@ -218,11 +223,14 @@ inline const std::vector& MCKinematicsReader::getTrackRefsBy { auto const& perEvent = mIndexedTrackRefs.at(source); if (perEvent.size() == 0) { - loadTrackRefsForSource(source); + initTrackRefsForSource(source); } if (static_cast(event) >= perEvent.size()) { reportMissingEvent("events of track references", source, event, perEvent.size()); } + if (!mTrackRefsLoaded[source][event]) { + loadTrackRefsForSourceAndEvent(source, event); + } return perEvent[event].getTruthArray(); } diff --git a/Steer/src/MCKinematicsReader.cxx b/Steer/src/MCKinematicsReader.cxx index 21024dba78368..42cd40c90c1ae 100644 --- a/Steer/src/MCKinematicsReader.cxx +++ b/Steer/src/MCKinematicsReader.cxx @@ -92,9 +92,15 @@ void MCKinematicsReader::loadTracksForSourceAndEvent(int source, int event) cons std::vector* loadtracks = nullptr; br->SetAddress(&loadtracks); br->GetEntry(event); - mTracks[source][event] = new std::vector; - *mTracks[source][event] = *loadtracks; - delete loadtracks; + // ROOT allocated the vector for us and we own it (we passed a pointer to nullptr): keep it instead of copying it + mTracks[source][event] = loadtracks; + br->ResetAddress(); // the branch must not refer to the stored vector (nor to the local pointer) any more + // free the decompressed baskets (~ the size of the event) if no later entry reads them, i.e. at the end of its cluster + auto clusterIt = br->GetTree()->GetClusterIterator(event); + clusterIt.Next(); + if (event + 1 >= clusterIt.GetNextEntry()) { + br->DropBaskets("all"); + } } } } @@ -105,6 +111,11 @@ void MCKinematicsReader::releaseTracksForSourceAndEvent(int source, int eventID) delete mTracks[source][eventID]; mTracks[source][eventID] = nullptr; } + // the track references of this event as well (reloaded on demand) + if (static_cast(eventID) < mTrackRefsLoaded.at(source).size() && mTrackRefsLoaded[source][eventID]) { + mIndexedTrackRefs[source][eventID] = o2::dataformats::MCTruthContainer(); + mTrackRefsLoaded[source][eventID] = false; + } } void MCKinematicsReader::loadHeadersForSource(int source) const @@ -129,31 +140,43 @@ void MCKinematicsReader::loadHeadersForSource(int source) const } } -void MCKinematicsReader::loadTrackRefsForSource(int source) const +void MCKinematicsReader::initTrackRefsForSource(int source) const { auto chain = mInputChains[source]; if (chain) { // todo: get name from NameConfig auto br = chain->GetBranch("TrackRefs"); if (br) { - std::vector* refs = nullptr; - br->SetAddress(&refs); mIndexedTrackRefs[source].resize(br->GetEntries()); - for (int event = 0; event < br->GetEntries(); ++event) { - br->GetEntry(event); - if (refs) { - // we convert the original flat vector into an indexed structure - initIndexedTrackRefs(*refs, mIndexedTrackRefs[source][event]); - delete refs; - refs = nullptr; - } - } + mTrackRefsLoaded[source].assign(br->GetEntries(), false); } else { LOG(warn) << "TrackRefs branch not found"; } } } +void MCKinematicsReader::loadTrackRefsForSourceAndEvent(int source, int event) const +{ + // todo: get name from NameConfig + auto br = mInputChains[source]->GetBranch("TrackRefs"); + std::vector* refs = nullptr; // allocated by ROOT, owned by us + br->SetAddress(&refs); + br->GetEntry(event); + if (refs) { + // we convert the original flat vector into an indexed structure + initIndexedTrackRefs(*refs, mIndexedTrackRefs[source][event]); + delete refs; + } + br->ResetAddress(); + // free the decompressed baskets if no later entry reads them, i.e. at the end of the cluster of this event + auto clusterIt = br->GetTree()->GetClusterIterator(event); + clusterIt.Next(); + if (event + 1 >= clusterIt.GetNextEntry()) { + br->DropBaskets("all"); + } + mTrackRefsLoaded[source][event] = true; +} + bool MCKinematicsReader::initFromDigitContext(o2::steer::DigitizationContext const* context) { if (mInitialized) { @@ -171,6 +194,7 @@ bool MCKinematicsReader::initFromDigitContext(o2::steer::DigitizationContext con mTracks.resize(mInputChains.size()); mHeaders.resize(mInputChains.size()); mIndexedTrackRefs.resize(mInputChains.size()); + mTrackRefsLoaded.resize(mInputChains.size()); // actual loading will be done only if someone asks // the first time for a particular source ... @@ -204,6 +228,7 @@ bool MCKinematicsReader::initFromKinematics(std::string_view name) mTracks.resize(1); mHeaders.resize(1); mIndexedTrackRefs.resize(1); + mTrackRefsLoaded.resize(1); mInitialized = true; return true;