BmnAlignTestTask.cxx

Go to the documentation of this file.

/*
    BM@N alignment routine
    BM@N experiment at NICA complex, JINR, 2025
 
    Department: Math & Soft Group of HEP lab
    Author:     Igor Polev, polev@jinr.int
 
    BmnAlignTestTask class implementation
*/
 
#include "BmnAlignTestTask.h"
 
#include "BmnAlignDefines.h"
#include "BmnHit.h"
#include "BmnRootIterator.h"
#include "BmnWeightedLine.h"
#include "fairlogger/Logger.h"
 
#include <TArrayI.h>
#include <TClonesArray.h>
#include <TDecompSVD.h>
#include <TFile.h>
#include <TMath.h>
#include <TMatrixD.h>
#include <TRandom.h>
#include <TStopwatch.h>
#include <TTree.h>
#include <TVector3.h>
#include <array>
#include <cmath>
#include <map>
#include <memory>
 
namespace
{
// truth solution
constexpr Double_t RANK_TOL{1.0E-10};        // singular values below this share of the largest are null
constexpr Double_t CONSTRAINT_TOL{1.0E-9};   // cm or rad, largest constraint violation the truth may keep
// vertex tracks
constexpr Int_t VERTEX_STEPS{3};   // Gauss-Newton steps of the vertex condition, see Make()
}   // anonymous namespace
 
/*
 *  The test task
 *  -------------
 *  Real tracks give the directions and positions, the truth solution gives the
 *  misalignment, the measurement model makes the hits. Make() goes in three
 *  stages, all on the data in RAM except the last one, which reads the source
 *  files once more to copy what the output needs.
 *
 *  1. Measurements on the real data, its hits taken in the base geometry (the
 *     aligner's last result, the nominal one without it):
 *     - the mean hit error of every element, reported and used as a noise option;
 *     - the vertex hit test of the aligner splits the tracks into vertex tracks
 *       (their vertex hit is accepted) and free tracks (rejected, or none);
 *     - z0, where the tracks meet the beam, see BmnAligner::FitTargetZ().
 *  2. The truth solution, see below.
 *  3. The hits. A free track keeps its free line: it misses the vertex as it
 *     does in the real data, and the aligner's vertex test has it to reject. A
 *     vertex track is the line closest to its real hits (in the metric of the
 *     free fit covariance) under one hard condition: the measurement model of
 *     the target, with its truth correction, reproduces the vertex hit plus the
 *     noise of the event. For Lin3 the condition is linear in the track and one
 *     step is exact; Lin6 adds terms quadratic in it, so VERTEX_STEPS
 *     Gauss-Newton steps are taken, the Jacobian from central differences, which
 *     are exact for a quadratic function. Every tracker hit is the measurement
 *     model at the source hit z plus Gauss noise in the element frame. The SiBT
 *     hits are copied, so the reader builds the same beam line and vertex hit.
 *
 *  z0
 *  --
 *  The truth puts the target at z0: its dZ gets z_c - z0 on top of the draw.
 *  The measurement model reads the vertex hit at z_c as the track at z_c - dZ
 *  (Lin3: X = track(z - dZ) + dX), so the tracks then meet the beam at z0.
 *
 *  The truth solution
 *  ------------------
 *  A test of the solver needs a truth the solver can reach. So the truth must
 *  satisfy every constraint the alignment will use, and be as likely as the
 *  truth sigmas - the aligner's presigmas - say. It is drawn from that prior,
 *  conditioned on the constraints.
 *
 *  The draw q is the own share of every parameter: q_k = Ref_k + Gaus(0, s_k),
 *  s_k = sigmaScale * Sigma_k. A parameter without a sigma is not drawn.
 *  The stored solution p holds effective leaf corrections, p = S q with
 *  S = [ I J ; 0 I ], and the constraints are written on p: C p = r, i.e.
 *  G q = r with G = C S.
 *
 *  1. A constraint with a single coefficient fixes that parameter, e.g. every
 *     rotation in a Lin3 alignment. It is not drawn but set to r / c.
 *  2. The rest is the conditional Gaussian draw, a projection in the metric of
 *     the prior covariance V = diag(s^2):
 *
 *         q' = q - K (G q - r),    K = V G^T (G V G^T)^+
 *
 *     Each parameter moves in proportion to its variance, so whatever was not
 *     drawn stays in place. The fixes of step 1 stay in G: on a leaf, they
 *     still have to cancel what the ancestors of the element bring.
 *
 *  Why q' is the draw under the constraints: q' is linear in q, so q' and G q
 *  are jointly Gaussian, and Cov(q', G q) = (I - K G) V G^T = 0. So they are
 *  independent: the law of q' does not depend on the value of G q, and where
 *  G q = r it is q itself. Hence law(q') = law(q | G q = r). q' is also the
 *  projection of q onto the constraints in the metric of V^-1, the minimum of
 *  (x - q)^T V^-1 (x - q) under G x = r. All constraint types are plain rows
 *  of C by now, so one step serves them all at once.
 */
 
Bool_t BmnAlignTestTask::Make(BmnAligner& aligner,
                              const char* dataPath,
                              const char* solutionPath,
                              Double_t sigmaScale,
                              ENoise noise,
                              const TVector2& stdDevXY)
{
    TStopwatch timer;
    LOG(state) << " > BmnAlignTestTask : preparing test task";
    if (!aligner.Initialized()) {
        LOG(error) << "<!> BmnAligner is not initialized";
        LOG(state) << " < BmnAlignTestTask : test task preparation terminated";
        return kFALSE;
    }
    if (!aligner.Hierarchy().Ready()) {
        LOG(error) << " <!> BmnAlignTestTask::Make: no parameter layout, PrepareData() must run first";
        return kFALSE;
    }
    if (sigmaScale < 0.0) {
        LOG(error) << " <!> BmnAlignTestTask::Make: sigmaScale " << sigmaScale << " is negative";
        return kFALSE;
    }
    if (kUserValues == noise && (stdDevXY.X() < 0.0 || stdDevXY.Y() < 0.0)) {
        LOG(error) << " <!> BmnAlignTestTask::Make: noise StdDev (" << stdDevXY.X() << ", " << stdDevXY.Y()
                   << ") is negative";
        return kFALSE;
    }
 
    const BmnAlignHierarchy& hier{aligner.Hierarchy()};
    const BmnDetectorModel& model{*aligner.DetectorModel()};
    const BmnMeasureModel& measure{*aligner.MeasureModel()};
    const Int_t nElem{hier.LeafCount()}, trgID{model.TargetID()}, trgIdx{trgID >= 0 ? model.Idx(trgID) : -1};
    const Bool_t withBeam{trgID >= 0};
    const Double_t zC{0.5 * (model.Target().zStart + model.Target().zEnd)};
 
    // prepare output file
    TFile outFile(dataPath, "RECREATE");
    if (outFile.IsZombie()) {
        LOG(error) << " <!> BmnAlignTestTask::Make: failed to create output file " << dataPath;
        return kFALSE;
    }
    outFile.cd();
 
    std::unique_ptr<BmnAlignResult> pResultCopy;
    auto RestoreResults = [&]() {
        if (pResultCopy)
            aligner.ResultLast() = *pResultCopy;
        else
            aligner.DropLastResult();   // it was null at the beginning of the method
    };
    VecSVectGL_t baseA(nElem);   // the base geometry, see the comment above
    if (aligner.Results().empty())
        aligner.AddBlankSolution();
    else {
        pResultCopy = std::make_unique<BmnAlignResult>(aligner.ResultLast());
        for (Int_t i = 0; i < nElem; i++)
            baseA[i] = aligner.ResultLast().A(i);
    }
 
    // Free lines through the tracker hits of the current track, in the base
    // geometry. Returns kTRUE if the track has a vertex hit, and gives its
    // position and weights (zero weights: disabled by the aligner's test).
    // Leaves the iterator past the last hit.
    auto FitFree = [&](BmnDataIterator& itr, BmnWeightedLine& lineX, BmnWeightedLine& lineY, SVect3& vertex,
                       SVect2& vertexW) {
        Bool_t hasVertex{kFALSE};
        for (; !itr.EndOfHits(); itr.NextHit()) {
            const Int_t detID{itr.HitDetectorID()};
            if (detID == trgID) {
                vertex = SVect3(itr.HitX(), itr.HitY(), itr.HitZ());
                vertexW = SVect2(itr.HitWx(), itr.HitWy());
                hasVertex = kTRUE;
                continue;
            }
            SVect3 hit{itr.HitX(), itr.HitY(), itr.HitZ()};
            measure.AlignHitApproximately(hit, baseA[model.Idx(detID)]);
            lineX.Add(hit(iZ), hit(iX), itr.HitWx());
            lineY.Add(hit(iZ), hit(iY), itr.HitWy());
        }
        return hasVertex;
    };
 
    // ==== stage 1: measurements on the real data ====
    BmnRamIterator& ram{*aligner.RamIterator()};
    VecDouble_t errX(nElem, 0.0), errY(nElem, 0.0);
    VecInt_t errN(nElem, 0);
    aligner.ProgressBar()->Init();
    for (ram.ResetAll(); !ram.EndOfTracks(); ram.NextTrack())
        for (; !ram.EndOfHits(); ram.NextHit()) {
            const Int_t idx{model.Idx(ram.HitDetectorID())};
            if (idx == trgIdx)
                continue;
            errX[idx] += 1.0 / std::sqrt(ram.HitWx());
            errY[idx] += 1.0 / std::sqrt(ram.HitWy());
            errN[idx]++;
        }
    aligner.ProgressBar()->Clear();
    LOG(info) << "    mean hit errors of the elements in the source data, x and y (cm):";
    for (Int_t i = 0; i < nElem; i++) {
        if (i == trgIdx)
            continue;
        errX[i] /= errN[i];   // PrepareData() refuses an element without hits
        errY[i] /= errN[i];
        LOG(info) << Form("      det %3d (%s) : %.2e  %.2e", i, model.DecodeID(model.EncodedID(i)).c_str(), errX[i],
                          errY[i]);
    }
 
    VecBool_t vertexTrack(aligner.TotalTracks(), kFALSE);
    Double_t z0{zC};
    if (withBeam) {
        // runs the vertex hit test too, whose hit states the flags below read
        Double_t z0Err;
        if (!aligner.FitTargetZ(z0, z0Err)) {
            z0 = zC;
            LOG(warning) << "  ! z0 is left at the target centre";
        }
        Int_t iTrack{0}, vertexTracks{0};
        aligner.ProgressBar()->Init();
        for (ram.ResetAll(); !ram.EndOfTracks(); ram.NextTrack(), iTrack++) {
            // the vertex hit comes first, and a disabled one has zero weights
            vertexTrack[iTrack] =
                !ram.EndOfHits() && ram.HitDetectorID() == trgID && ram.HitWx() > 0.0 && ram.HitWy() > 0.0;
            vertexTracks += vertexTrack[iTrack];
        }
        aligner.ProgressBar()->Clear();
        LOG(info) << Form("    %d vertex tracks, %zu free tracks (vertex hit rejected or absent)", vertexTracks,
                          vertexTrack.size() - vertexTracks);
    }   // if(withBeam)
 
    // ==== stage 2: the truth solution, see the comment above ====
    const BmnAlignPriors& priors{aligner.Priors()};
    const Int_t nFull{hier.AllParams()}, nCnstr{priors.ConstCount()};
    TMatrixD mS(nFull, nFull);
    mS.UnitMatrix();
    if (hier.CompositeParams())   // SetSub refuses the empty J of a flat model
        mS.SetSub(0, hier.CompositeParamBase(), hier.J());
 
    TVectorD vQ(nFull),   // own shares
        vVar(nFull);      // their prior variance, zero where nothing is drawn
    for (const ParamRef& p : hier.Layout().All()) {
        const Double_t sigma{sigmaScale * priors.SigmaAt(p.full)};
        if (sigma <= 0.0)
            continue;
        vQ(p.full) = priors.ref[p.full] + gRandom->Gaus(0.0, sigma);
        vVar(p.full) = sigma * sigma;
    }
    if (withBeam)
        vQ(hier.Layout().ByNode()[trgIdx].params[idZ].full) += zC - z0;
 
    // RMS of the effective element corrections and of the composite own shares,
    // per correction. The target is left out: its dZ carries z_c - z0.
    auto LogRms = [&](const TVectorD& q, const char* stage) {
        const TVectorD p{mS * q};
        SVectGL leaf, comp;
        for (const ParamRef& r : hier.Layout().All()) {
            if (!r.IsComposite() && r.leaf == trgIdx)
                continue;   // works for targetless model also
            (r.IsComposite() ? comp : leaf)(r.prm) += p(r.full) * p(r.full);
        }
        const Int_t nLeaf{nElem - withBeam}, nComp{hier.CompositeCount()};
        auto Row = [](const SVectGL& sum, Int_t n) {
            TString row;
            for (Int_t k = 0; k < BMN_GLOBAL_PARAMS_PD; k++)
                row += Form(" %10.3e", n ? std::sqrt(sum(k) / n) : 0.0);
            return row;
        };
        LOG(info) << Form("      %-12s elements  ", stage) << Row(leaf, nLeaf);
        LOG_IF(info, nComp) << Form("      %-12s composites", stage) << Row(comp, nComp);
    };
    TString head;
    for (Int_t k = 0; k < BMN_GLOBAL_PARAMS_PD; k++)
        head += Form(" %10s", BmnParamLayout::ParamName(k));
    LOG(info) << Form("      %-12s %-10s", "truth RMS", "(cm, rad)") << head;
    LogRms(vQ, "drawn");
 
    VecBool_t isFixed(nFull, kFALSE);
    for (Int_t j = 0; j < nCnstr; j++) {
        const TMatrixDRow_const c(priors.constraints, j);
        Int_t k{-1}, terms{0};
        for (Int_t i = 0; i < nFull; i++)
            if (0.0 != c(i)) {
                k = i;
                terms++;
            }
        if (1 != terms)
            continue;
        vQ(k) = priors.rhs(j) / c(k);
        vVar(k) = 0.0;
        isFixed[k] = kTRUE;
    }
    Int_t drawn{0}, fixed{0};
    for (Int_t i = 0; i < nFull; i++) {
        if (vVar(i) > 0.0)
            drawn++;
        if (isFixed[i])
            fixed++;
    }
 
    if (nCnstr) {
        const TVectorD& vR{priors.rhs};
        const TMatrixD mG{priors.constraints, TMatrixD::kMult, mS};
        TMatrixD mGV{mG};                                       // G V
        mGV.NormByRow(vVar, "M");                               // despite the name: column i times vVar(i)
        const TMatrixD mM{mGV, TMatrixD::kMultTranspose, mG};   // G V G^T
        // all-zero when nothing was drawn, and SVD does not survive that
        if (mM.NonZeros()) {
            TVectorD vRes{mG * vQ};
            vRes -= vR;
            TDecompSVD svd(mM);
            if (!svd.Decompose()) {
                LOG(error) << " <!> BmnAlignTestTask::Make: SVD failed, the truth cannot be projected"
                              " onto the constraints";
                RestoreResults();
                return kFALSE;
            }
            // pseudo-inverse: redundant constraints leave null singular values
            const TVectorD& sig{svd.GetSig()};   // TDecompSVD returns them in decreasing order
            TVectorD vX{TMatrixD(TMatrixD::kTransposed, svd.GetU()) * vRes};
            for (Int_t j = 0; j < nCnstr; j++)
                vX(j) = sig(j) > sig(0) * RANK_TOL ? vX(j) / sig(j) : 0.0;
            vX = svd.GetV() * vX;
            vQ -= TMatrixD(TMatrixD::kTransposed, mGV) * vX;
        }
        TVectorD vViol{mG * vQ};
        vViol -= vR;
        Int_t violated{0};
        for (Int_t j = 0; j < nCnstr; j++)
            if (std::fabs(vViol(j)) > CONSTRAINT_TOL) {
                LOG(error) << " <!> BmnAlignTestTask::Make: constraint #" << j + 1 << " cannot hold, violated by "
                           << vViol(j) << " - it rests on parameters without a truth sigma";
                violated++;
            }
        if (violated) {
            RestoreResults();
            return kFALSE;
        }
    }
    LogRms(vQ, "constrained");
 
    BmnAlignResult& truth{aligner.ResultLast()};
    truth.A() = mS * vQ;
    truth.SeparateA();
    LOG(info) << "    truth drawn: " << drawn << " parameter(s) from the truth sigmas times " << sigmaScale << ", "
              << fixed << " fixed by single-parameter constraints, all " << nCnstr << " constraint(s) hold";
    LOG_IF(warning, !drawn && sigmaScale > 0.0)
        << "  ! no parameter carries a truth sigma: the truth is the nominal geometry";
    if (withBeam) {
        const SVectGL& aTrg{truth.A(trgIdx)};
        LOG(info) << Form("    target truth : dX %+.3e, dY %+.3e, dZ %+.3e cm", aTrg(idX), aTrg(idY), aTrg(idZ));
    }
 
    // ==== stage 3: the synthetic hits ====
    const Bool_t noiseX{kUserValues != noise || 0.0 != stdDevXY.X()},
        noiseY{kUserValues != noise || 0.0 != stdDevXY.Y()};
    TClonesArray trackData("CbmStsTrack", 200), hitData("CbmStsHit", 200 * model.Alignables()), beamData("BmnHit"),
        *tdataPtr{&trackData}, *hdataPtr{&hitData}, *bdataPtr{&beamData};
    BmnRootIterator& root{*aligner.RootIterator()};
    // Note ------------------------------------------------------------------------
    // When copying a TTree structure with CloneTree(0) if data is governed by
    // TChain holding more than one file, modification of data will be lost and
    // replaced by original data.
    // Cause -----------------------------------------------------------------------
    // CloneTree(0) registers the clone for TChain::Notify() callbacks.
    // When the chain transitions between files (including on ResetAll/Restart),
    // Notify() resets the clone's branch addresses back to the chain's own buffers,
    // silently undoing SetBranchAddress overrides. This causes Fill() to write the
    // original (unmodified) chain data instead of our modified trackData/hitData.
    // Conclusion ------------------------------------------------------------------
    // Can't use CloneTree(0) to copy structure of TTree from original data,
    // must create an independent TTree instead.
    // Aknowledgement --------------------------------------------------------------
    // Thanks to Sonnet 4.6 for this non-trivial issue diagnostic.
    TTree outTree(root.GetChain().GetTree()->GetName(), root.GetChain().GetTree()->GetTitle());
    outTree.SetDirectory(&outFile);
    outTree.Branch(BMN_TRACKS_BRANCH_NAME, "TClonesArray", &tdataPtr);
    outTree.Branch(BMN_HITS_BRANCH_NAME, "TClonesArray", &hdataPtr);
    if (withBeam)
        outTree.Branch(BMN_BEAM_BRANCH_NAME, "TClonesArray", &bdataPtr);
 
    std::map<Int_t, std::array<Double_t, 3>> sibtErr;   // per SiBT station: sum of dx, of dy, hit count
    Int_t iTrack{0}, freeWithVertex{0};
    Double_t vertexNoiseX{0.0},   // beam error of the current event, common to its tracks
        vertexNoiseY{0.0};
    aligner.ProgressBar()->Init();
    root.ResetAll();
    while (!root.EndOfTracks()) {   // `while` instead of `for` is necessary here, see comment below
        if (withBeam && root.IsNewEvent()) {
            const TClonesArray* beamHits{root.GetBeamHits()};
            for (Int_t i = 0; i < beamHits->GetEntriesFast(); i++) {
                const BmnHit& hit{*static_cast<const BmnHit*>(beamHits->At(i))};
                new (beamData[i]) BmnHit(hit);
                std::array<Double_t, 3>& err{sibtErr[hit.GetStation()]};
                err[0] += hit.GetDx();
                err[1] += hit.GetDy();
                err[2] += 1.0;
            }
            vertexNoiseX = noiseX ? gRandom->Gaus(0.0, std::sqrt(root.VertexVarX())) : 0.0;
            vertexNoiseY = noiseY ? gRandom->Gaus(0.0, std::sqrt(root.VertexVarY())) : 0.0;
        }
 
        // the true track
        BmnWeightedLine lineX, lineY;
        SVect3 vertex;
        SVect2 vertexW;
        const Bool_t hasVertex{FitFree(root, lineX, lineY, vertex, vertexW)};
        if (lineX.Det() <= 0.0 || lineY.Det() <= 0.0) {
            LOG(error) << " <!> BmnAlignTestTask::Make: track #" << iTrack << " has too few hits for a straight line";
            RestoreResults();
            return kFALSE;
        }
        const SVectLC alphaFree{lineX.At0(), lineX.Slope(), lineY.At0(), lineY.Slope()};
        SVectLC alpha{alphaFree};
        if (vertexTrack[iTrack]) {
            const SVectGL& aTrg{truth.A(trgIdx)};
            const Double_t zV{vertex(iZ)};
            auto Vertex = [&](const SVectLC& a) { return SVect2(measure.X(zV, aTrg, a), measure.Y(zV, aTrg, a)); };
            const SVect2 vB{vertex(iX) + vertexNoiseX, vertex(iY) + vertexNoiseY};
            SMatrLC mCov;   // of the free lines, x and y apart
            mCov(iX0, iX0) = lineX.Var(0.0);
            mCov(iX0, iX1) = lineX.Cov();   // symmetric storage: (iX1, iX0) is the same element
            mCov(iX1, iX1) = lineX.SlopeVar();
            mCov(iY0, iY0) = lineY.Var(0.0);
            mCov(iY0, iY1) = lineY.Cov();
            mCov(iY1, iY1) = lineY.SlopeVar();
            for (Int_t step = 0; step < VERTEX_STEPS; step++) {
                RM::SMatrix<Double_t, 2, BMN_LOCAL_PARAMS_PT> mA;   // Jacobian of the condition at alpha
                for (Int_t j = 0; j < BMN_LOCAL_PARAMS_PT; j++) {
                    SVectLC up{alpha}, down{alpha};
                    up(j) += 1.0;
                    down(j) -= 1.0;
                    const SVect2 d{0.5 * (Vertex(up) - Vertex(down))};
                    mA(0, j) = d(0);
                    mA(1, j) = d(1);
                }
                const RM::SMatrix<Double_t, BMN_LOCAL_PARAMS_PT, 2> mCAt{mCov * RM::Transpose(mA)};
                RM::SMatrix<Double_t, 2, 2, RM::MatRepSym<Double_t, 2>> mM{RM::Similarity(mA, mCov)};   // A C A^T
                if (!mM.InvertChol()) {
                    LOG(error) << " <!> BmnAlignTestTask::Make: track #" << iTrack
                               << " cannot be put through the vertex";
                    RestoreResults();
                    return kFALSE;
                }
                // the free line projected onto the condition linearized at alpha
                alpha = alphaFree - mCAt * (mM * (mA * (alphaFree - alpha) + Vertex(alpha) - vB));
            }
        } else if (hasVertex)
            freeWithVertex++;
 
        // the synthetic hits
        Int_t arrIdx{0};
        TArrayI hitArray(root.HitsInTrack() - root.HasVertexHit());
        for (root.ResetHits(); !root.EndOfHits(); root.NextHit()) {
            const BmnRootIterator::Hit_t* srcHit{root.GetHit()};
            if (!srcHit)
                continue;   // the vertex hit: the reader builds it again from the SiBT hits
            const Int_t hitIdx{hitData.GetEntriesFast()}, hitDet{srcHit->GetDetectorID()},
                detIdx{model.IdxFromHit(hitDet)};
            const Double_t z0hit{srcHit->GetZ()},
                dx{kHitErrors == noise      ? srcHit->GetDx()
                   : kElementMeans == noise ? errX[detIdx]
                                            : stdDevXY.X()},
                dy{kHitErrors == noise      ? srcHit->GetDy()
                   : kElementMeans == noise ? errY[detIdx]
                                            : stdDevXY.Y()};
            TVector3 pos{// not const: CbmStsHit takes them by non-const reference
                         measure.X(z0hit, truth.A(detIdx), alpha) + (dx ? gRandom->Gaus(0, dx) : 0.0),
                         measure.Y(z0hit, truth.A(detIdx), alpha) + (dy ? gRandom->Gaus(0, dy) : 0.0), z0hit},
                dpos{dx ? dx : BmnHitRecord::MIN_HIT_ERROR, dy ? dy : BmnHitRecord::MIN_HIT_ERROR, 0.0};
            new (hitData[hitIdx]) CbmStsHit(hitDet, pos, dpos, 0.0, 0, 0);
            hitArray[arrIdx++] = hitIdx;
        }
 
        // the synthetic track
        BmnRootIterator::Track_t* srcTrack{const_cast<BmnRootIterator::Track_t*>(root.GetTrack())};
        CbmStsTrack* newTrack = new (trackData[trackData.GetEntriesFast()]) CbmStsTrack();
        newTrack->SetStsHits(hitArray);
        newTrack->SetTrkID(srcTrack->GetTrkID());
        newTrack->SetPidHypo(srcTrack->GetPidHypo());
        newTrack->SetParamFirst(*srcTrack->GetParamFirst());
        newTrack->SetParamLast(*srcTrack->GetParamLast());
        newTrack->SetFlag(srcTrack->GetFlag());
        newTrack->SetChi2(srcTrack->GetChi2());
        newTrack->SetNDF(srcTrack->GetNDF());
        newTrack->SetB(srcTrack->GetB());
        newTrack->SetStsEv(srcTrack->GetStsEv());
 
        root.NextTrack();
        iTrack++;
        // flush accumulated event if the next track opens a new one or there is none.
        // CAUTION: the test looks one track ahead - NextTrack() must stay right before it
        if (root.IsNewEvent() || root.EndOfTracks()) {
            outFile.cd();   // just in case...
            outTree.Fill();
            trackData.Delete();   // not using Clear() to avoid memory leaks:
            hitData.Delete();     // maybe it is safe now, but can't be sure in future changes
            beamData.Delete();
        }
    }   // loop over tracks
    aligner.ProgressBar()->Clear();
 
    for (const auto& [station, err] : sibtErr)
        LOG(info) << Form("    SiBT station %d : mean hit errors x %.2e, y %.2e cm", station, err[0] / err[2],
                          err[1] / err[2]);
    LOG_IF(info, withBeam) << "    " << freeWithVertex
                           << " free tracks take a vertex hit in the task:"
                              " the aligner's vertex test should reject about as many";
 
    outFile.Write("", TObject::kOverwrite);
    outFile.Close();
    LOG(info) << "    test data successfully written to " << dataPath;
    if (aligner.SaveSolution(solutionPath)) {
        LOG(info) << "    test solution successfully written to  " << solutionPath;
    } else {
        LOG(error) << " <!> BmnAlignTestTask::Make: failed to save solution to " << solutionPath;
        RestoreResults();
        return kFALSE;
    }
 
    RestoreResults();
    LOG(state) << Form(" < BmnAlignTestTask : test task prepared in %.1fs", timer.RealTime());
    return kTRUE;
}