diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h index 80a47c5ce6405..fe33e3baa5c7c 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertexExt.h @@ -26,10 +26,17 @@ struct PrimaryVertexExt : public PrimaryVertex { using PrimaryVertex::PrimaryVertex; std::array nSrc{}; // N contributors for each source type int VtxID = -1; // original vtx ID - float FT0Amp = -1; // amplitude of closest FT0 trigger - float FT0A = -1; // amplitude of the A side + float FT0A = -1; // amplitude of closest FT0 A side + float FT0C = -1; // amplitude of closest FT0 C side float FT0Time = -1.; // time of closest FT0 trigger - + float rmsT = 0; + float rmsZ = 0; + float rmsTW = 0; + float rmsZW = 0; + float rmsT0 = 0; // w/o ITS + float rmsTW0 = 0; // w/o ITS + float tMAD = 0; + float zMAD = 0; int getNSrc(int i) const { return nSrc[i]; } #ifndef GPUCA_ALIGPUCODE @@ -37,7 +44,7 @@ struct PrimaryVertexExt : public PrimaryVertex { std::string asString() const; #endif - ClassDefNV(PrimaryVertexExt, 2); + ClassDefNV(PrimaryVertexExt, 4); }; #ifndef GPUCA_ALIGPUCODE diff --git a/DataFormats/Reconstruction/src/PrimaryVertexExt.cxx b/DataFormats/Reconstruction/src/PrimaryVertexExt.cxx index 3d8ec6c5fa7c2..6065f04a3bc1a 100644 --- a/DataFormats/Reconstruction/src/PrimaryVertexExt.cxx +++ b/DataFormats/Reconstruction/src/PrimaryVertexExt.cxx @@ -24,7 +24,7 @@ using GTrackID = o2::dataformats::GlobalTrackID; std::string PrimaryVertexExt::asString() const { - auto str = o2::utils::Str::concat_string(PrimaryVertex::asString(), fmt::format("VtxID={} FT0A/C={}/{} FT0T={}", VtxID, FT0A, FT0Amp - FT0A, FT0Time)); + auto str = o2::utils::Str::concat_string(PrimaryVertex::asString(), fmt::format("VtxID={} FT0A/C={}/{} FT0T={}", VtxID, FT0A, FT0C, FT0Time)); for (int i = 0; i < GTrackID::Source::NSources; i++) { if (getNSrc(i) > 0) { str += fmt::format(" {}={}", GTrackID::getSourceName(i), getNSrc(i)); diff --git a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt index bfa253dd5980a..21142374ba25f 100644 --- a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt @@ -9,6 +9,8 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +# add_compile_options(-O0 -g -fPIC) + o2_add_library(GlobalTrackingStudy SOURCES src/TPCTrackStudy.cxx src/TrackingStudy.cxx diff --git a/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx b/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx index ce2d53d71c385..94b9dd311507a 100644 --- a/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx +++ b/Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx @@ -37,6 +37,7 @@ #include "ReconstructionDataFormats/VtxTrackRef.h" #include "ReconstructionDataFormats/DCA.h" #include "Steer/MCKinematicsReader.h" +#include "MathUtils/fit.h" namespace o2::trackstudy { @@ -74,7 +75,7 @@ class TrackingStudySpec : public Task std::unique_ptr mDBGOutVtx; float mITSROFrameLengthMUS = 0.; int mMaxNeighbours = 3; - float mMaxVTTimeDiff = 25.; // \mus + float mMaxVTTimeDiff = 80.; // \mus GTrackID::mask_t mTracksSrc{}; o2::steer::MCKinematicsReader mcReader; // reader of MC information }; @@ -125,14 +126,19 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData) int nv = vtxRefs.size() - 1; o2::dataformats::PrimaryVertexExt pveDummy; std::vector pveVec; + float tBiasITS = o2::itsmft::DPLAlpideParam::Instance().roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS; + std::vector dtvec, dzvec; for (int iv = 0; iv < nv; iv++) { LOGP(debug, "processing PV {} of {}", iv, nv); - const auto& vtref = vtxRefs[iv]; auto pv = pvvec[iv]; auto& pve = pveVec.emplace_back(); static_cast(pve) = pv; + dtvec.clear(); + dzvec.clear(); + dtvec.reserve(pv.getNContributors()); + dzvec.reserve(pv.getNContributors()); float bestTimeDiff = 1000, bestTime = -999; int bestFTID = -1; @@ -152,10 +158,15 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData) } if (bestFTID >= 0) { pve.FT0A = FITInfo[bestFTID].getTrigger().getAmplA(); - pve.FT0Amp = pve.FT0A + FITInfo[bestFTID].getTrigger().getAmplC(); + pve.FT0C = FITInfo[bestFTID].getTrigger().getAmplC(); pve.FT0Time = FITInfo[bestFTID].getInteractionRecord().differenceInBCMUS(recoData.startIR); } pve.VtxID = iv; + float meanT = 0, meanZ = 0, rmsT = 0, rmsZ = 0; + float meanTW = 0, meanZW = 0, rmsTW = 0, rmsZW = 0, WT = 0, WZ = 0; + float meanT0 = 0, rmsT0 = 0; + float meanTW0 = 0, rmsTW0 = 0, WT0 = 0; + int nContAdd = 0, nContAdd0 = 0, ntITS = 0; for (int is = 0; is < GTrackID::NSources; is++) { bool skipTracks = (!GTrackID::getSourceDetectorsMask(is)[GTrackID::ITS] || !mTracksSrc[is] || !recoData.isTrackSourceLoaded(is)); int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is); @@ -177,13 +188,97 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData) } float ttime = 0, ttimeE = 0; recoData.getTrackTime(vid, ttime, ttimeE); + bool acceptForPV0 = pvCont; + if (vid.getSource() == GTrackID::ITS) { + ttimeE *= mITSROFrameLengthMUS; + ttime += ttimeE + tBiasITS; + ttimeE *= 1. / sqrt(3); + if (++ntITS > 0) { // do not allow ITS in the MAD + acceptForPV0 = false; + } + } else if (vid.getSource() == GTrackID::TPC) { + ttimeE *= o2::constants::lhc::LHCBunchSpacingMUS * 8; + } + if (pvCont) { + float dt = ttime - pve.getTimeStamp().getTimeStamp(); + float tW = 1. / (ttimeE * ttimeE), zW = 1. / trc.getSigmaZ2(); + dzvec.push_back(dca.getZ()); + WT += tW; + WZ += zW; + meanT += dt; + meanTW += dt * tW; + meanZ += dca.getZ(); + meanZW += dca.getZ() * zW; + + rmsT += dt * dt; + rmsTW += dt * dt * tW; + rmsZ += dca.getZ() * dca.getZ(); + rmsZW += dca.getZ() * dca.getZ() * zW; + nContAdd++; + if (acceptForPV0) { + dtvec.push_back(dt); + WT0 += tW; + meanT0 += dt; + meanTW0 += dt * tW; + rmsT0 += dt * dt; + rmsTW0 += dt * dt * tW; + nContAdd0++; + } + LOGP(debug, "dt={} dz={}, tW={}, zW={} t={} tE={} {}", dt, dca.getZ(), tW, zW, ttime, ttimeE, vid.asString()); + } (*mDBGOut) << "dca" << "tfID=" << TFCount << "ttime=" << ttime << "ttimeE=" << ttimeE << "gid=" << vid << "pv=" << pv << "trc=" << trc << "pvCont=" << pvCont << "ambig=" << ambig << "dca=" << dca << "xmin=" << xmin << "\n"; } } + + if (nContAdd) { + rmsT /= nContAdd; + rmsZ /= nContAdd; + meanT /= nContAdd; + meanZ /= nContAdd; + pve.rmsT = (rmsT - meanT * meanT); + pve.rmsT = pve.rmsT > 0 ? std::sqrt(pve.rmsT) : 0; + pve.rmsZ = rmsZ - meanZ * meanZ; + pve.rmsZ = pve.rmsZ > 0 ? std::sqrt(pve.rmsZ) : 0; + } + if (nContAdd0) { + rmsT0 /= nContAdd0; + meanT0 /= nContAdd0; + pve.rmsT0 = (rmsT0 - meanT0 * meanT0); + pve.rmsT0 = pve.rmsT0 > 0 ? std::sqrt(pve.rmsT0) : 0; + } + if (WT0 > 0) { + rmsTW0 /= WT0; + meanTW0 /= WT0; + pve.rmsTW0 = (rmsTW0 - meanTW0 * meanTW0); + pve.rmsTW0 = pve.rmsTW0 > 0 ? std::sqrt(pve.rmsTW0) : 0; + } + // + if (WT > 0 && WZ > 0) { + rmsTW /= WT; + meanTW /= WT; + pve.rmsTW = (rmsTW - meanTW * meanTW); + pve.rmsTW = pve.rmsTW > 0 ? std::sqrt(pve.rmsTW) : 0; + rmsZW /= WZ; + meanZW /= WZ; + pve.rmsZW = rmsZW - meanZW * meanZW; + pve.rmsZW = pve.rmsZ > 0 ? std::sqrt(pve.rmsZ) : 0; + } + pve.tMAD = o2::math_utils::MAD2Sigma(dtvec.size(), dtvec.data()); + pve.zMAD = o2::math_utils::MAD2Sigma(dzvec.size(), dzvec.data()); } int nvtot = mMaxNeighbours < 0 ? -1 : (int)pveVec.size(); + + auto insSlot = [maxSlots = mMaxNeighbours](std::vector& vc, float v, int slot, std::vector& vid, int id) { + for (int i = maxSlots - 1; i > slot; i--) { + std::swap(vc[i], vc[i - 1]); + std::swap(vid[i], vid[i - 1]); + } + vc[slot] = v; + vid[slot] = id; + }; + for (int cnt = 0; cnt < nvtot; cnt++) { const auto& pve = pveVec[cnt]; float tv = pve.getTimeStamp().getTimeStamp(); @@ -205,16 +300,14 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData) } for (int i = 0; i < mMaxNeighbours; i++) { if (dT[i] > dtime) { - dT[i] = dtime; - idT[i] = cntM; + insSlot(dT, dtime, i, idT, cntM); break; } } auto dz = std::abs(pve.getZ() - vt.getZ()); for (int i = 0; i < mMaxNeighbours; i++) { if (dZ[i] > dz) { - dZ[i] = dz; - idZ[i] = cntM; + insSlot(dZ, dz, i, idZ, cntM); break; } } @@ -227,16 +320,14 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData) } for (int i = 0; i < mMaxNeighbours; i++) { if (dT[i] > dtime) { - dT[i] = dtime; - idT[i] = cntP; + insSlot(dT, dtime, i, idT, cntP); break; } } auto dz = std::abs(pve.getZ() - vt.getZ()); for (int i = 0; i < mMaxNeighbours; i++) { if (dZ[i] > dz) { - dZ[i] = dz; - idZ[i] = cntP; + insSlot(dZ, dz, i, idZ, cntP); break; } } @@ -302,8 +393,8 @@ DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mas outputs, AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, srcTracks, useMC)}, Options{ - {"max-vtx-neighbours", VariantType::Int, 2, {"Max PV neighbours fill, no PV study if < 0"}}, - {"max-vtx-timediff", VariantType::Float, 25.f, {"Max PV time difference to consider"}}, + {"max-vtx-neighbours", VariantType::Int, 3, {"Max PV neighbours fill, no PV study if < 0"}}, + {"max-vtx-timediff", VariantType::Float, 90.f, {"Max PV time difference to consider"}}, }}; }