vertex.C

Go to the documentation of this file.

/*
    BM@N alignment routine
    BM@N experiment at NICA complex, JINR, 2026
 
    Department: Math & Soft Group of HEP lab
    Author:     Claude (Opus 5)
    Supervisor: Igor Polev, polev@jinr.int
 
    Primary vertex study, phase 0 of the TRG expansion (see Expansion.md):
    how well the beam line measured by SiBT predicts where the STS tracks start,
    measured before any of it is built into the aligner
*/
 
// headers below are for IDE syntax to work and they will be
// turned off during execution in ROOT macro mode
#ifndef __CLING__
#include "BmnAlignModels.h"   // IWYU pragma: keep
#include "BmnDataFilters.h"
#include "BmnDataReader.h"
#include "BmnHit.h"
#include "BmnHitRecord.h"
#include "BmnWeightedLine.h"
#include "CbmStsHit.h"
#include "CbmStsTrack.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TFile.h"
#include "TMath.h"
#include "TProfile.h"
#include "TTree.h"
#include "TTreeReaderValue.h"
 
#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>
#endif
 
#include "fairlogger/Logger.h"
 
/**
 * @brief Compare the SiBT beam line with the STS tracks at the target centre.
 *
 * For every event with one hit in each SiBT station the beam is fitted as a
 * straight line and taken at the target centre z_c. Every STS track the aligner
 * would accept is fitted the way the aligner fits it - a weighted straight line
 * per projection - and taken at the same z_c. The per-event ntuple is written
 * out for interactive study, and a summary answers the phase 0 questions: the
 * real beam error, where the tracks meet the beam and how far they scatter around
 * it, the window that selects vertex tracks, and how many tracks of one event
 * share its vertex hit (see Expansion.md).
 *
 * The residual is r = beam - track at z_c. A track from z_v = z_c + dz gives
 * r = (track slope - beam slope) * dz. A z offset of the tracker relative to the
 * target has the same form, so the dz fitted from the slope dependence is the
 * vertex z in the tracker frame, not the interaction depth. Before the alignment
 * the station halves are shifted against each other, so tracks going to either
 * side of the beam get a line of their own.
 *
 * @param fileListPath text file listing the DST files, '#' lines are ignored
 * @param outPath      .root file for the ntuple and the profiles
 * @param maxEvents    number of events to read, 0 for all of them
 *
 * @return 0 on success, -1 on a setup failure
 */
int vertex(const char* fileListPath = "/home/igor/bmnroot/align_new/cfg/run9/9555.noa_rlerr.sample.list",
           const char* outPath = "vertex.root",
           Long64_t maxEvents = 1000)
{
    // the collision target installed in the run, as a FairRoot ASCII geometry file
    constexpr const char* TARGET_GEO{"/home/igor/bmnroot/geometry/target_Sn.geo"};
    constexpr Int_t SIBT_STATIONS{3}, SIBT_MIDDLE{1};   // stations are numbered along the beam
    constexpr Double_t WINDOW_SIGMAS{5.0},
        // this close to the beam direction in both projections a track
        // does not point to the vertex: its residuals spread by centimetres
        BEAM_LIKE_SLOPE{0.02};
 
    fair::Logger::SetConsoleColor(true);
    fair::Logger::SetConsoleSeverity("info");
 
    BmnRun9DM detectorModel;
    if (!detectorModel.Init(TARGET_GEO))
        return -1;
    const Double_t zStart{detectorModel.Target().zStart}, zEnd{detectorModel.Target().zEnd}, zc{0.5 * (zStart + zEnd)},
        halfT{0.5 * (zEnd - zStart)};
 
    BmnDataReader data{fileListPath};
    TTreeReader& reader{data.Reader()};
    TTreeReaderValue<TClonesArray> sibtHits{reader, "BmnSiBTHit"}, stsTracks{reader, "StsVector"},
        stsHits{reader, "StsHit"};
 
    TFile out{outPath, "RECREATE"};
    if (out.IsZombie())
        return -1;
    // owned by the file, which deletes it on Close()
    TTree* tree{new TTree("vertex", "SiBT beam line vs straight STS tracks at the target centre")};
    Float_t beamX, beamY, beamErrX, beamErrY, beamSlopeX, beamSlopeY, beamPullX, beamPullY;
    Int_t tracksAll;
    std::vector<Float_t> slopeX, slopeY, resX, resY, errX, errY;
    tree->Branch("beamX", &beamX);
    tree->Branch("beamY", &beamY);
    tree->Branch("beamErrX", &beamErrX);
    tree->Branch("beamErrY", &beamErrY);
    tree->Branch("beamSlopeX", &beamSlopeX);
    tree->Branch("beamSlopeY", &beamSlopeY);
    tree->Branch("beamPullX", &beamPullX);
    tree->Branch("beamPullY", &beamPullY);
    tree->Branch("tracksAll", &tracksAll);
    tree->Branch("slopeX", &slopeX);
    tree->Branch("slopeY", &slopeY);
    tree->Branch("resX", &resX);
    tree->Branch("resY", &resY);
    tree->Branch("errX", &errX);
    tree->Branch("errY", &errY);
 
    Long64_t events{0}, eventsWithBeam{0}, tracksKept{0}, ndf{0};
    Double_t chi2X{0.0}, chi2Y{0.0};
    while ((maxEvents <= 0 || events < maxEvents) && reader.Next()) {
        events++;
        if (SIBT_STATIONS != sibtHits->GetEntriesFast())
            continue;
 
        BmnWeightedLine bx, by;
        UInt_t stations{0};
        Double_t midZ{0.0}, midX{0.0}, midY{0.0};
        for (Int_t i{0}; i < SIBT_STATIONS; i++) {
            const BmnHit& hit{*static_cast<BmnHit*>(sibtHits->At(i))};
            if (hit.GetStation() >= 0 && hit.GetStation() < SIBT_STATIONS)
                stations |= 1u << hit.GetStation();
            if (SIBT_MIDDLE == hit.GetStation()) {
                midZ = hit.GetZ();
                midX = hit.GetX();
                midY = hit.GetY();
            }
            bx.Add(hit.GetZ(), hit.GetX(), BmnHitRecord::Weight(hit.GetDx()));
            by.Add(hit.GetZ(), hit.GetY(), BmnHitRecord::Weight(hit.GetDy()));
        }
        if (stations != (1u << SIBT_STATIONS) - 1)
            continue;
        eventsWithBeam++;
 
        beamX = bx.At(zc);
        beamY = by.At(zc);
        beamErrX = std::sqrt(bx.Var(zc));
        beamErrY = std::sqrt(by.Var(zc));
        beamSlopeX = bx.Slope();
        beamSlopeY = by.Slope();
        // three points leave one degree of freedom: all residuals are proportional
        // to a single pull, chi2 is its square, and the middle hit gives its sign
        beamPullX = std::copysign(std::sqrt(std::max(bx.Chi2(), 0.0)), midX - bx.At(midZ));
        beamPullY = std::copysign(std::sqrt(std::max(by.Chi2(), 0.0)), midY - by.At(midZ));
        tracksAll = stsTracks->GetEntriesFast();
        slopeX.clear();
        slopeY.clear();
        resX.clear();
        resY.clear();
        errX.clear();
        errY.clear();
 
        for (Int_t i{0}; i < tracksAll; i++) {
            const CbmStsTrack& track{*static_cast<CbmStsTrack*>(stsTracks->At(i))};
            BmnWeightedLine tx, ty;
            for (Int_t h{0}; h < track.GetNStsHits(); h++) {
                const CbmStsHit& hit{*static_cast<CbmStsHit*>(stsHits->At(track.GetStsHitIndex(h)))};
                if (!BmnAcceptHit(hit, detectorModel))
                    continue;
                tx.Add(hit.GetZ(), hit.GetX(), BmnHitRecord::Weight(hit.GetDx()));
                ty.Add(hit.GetZ(), hit.GetY(), BmnHitRecord::Weight(hit.GetDy()));
            }
            if (!BmnAcceptTrack(track, tx.N()))
                continue;
 
            chi2X += tx.Chi2();
            chi2Y += ty.Chi2();
            ndf += tx.N() - 2;
            slopeX.push_back(tx.Slope());
            slopeY.push_back(ty.Slope());
            resX.push_back(bx.At(zc) - tx.At(zc));
            resY.push_back(by.At(zc) - ty.At(zc));
            errX.push_back(std::sqrt(tx.Var(zc)));
            errY.push_back(std::sqrt(ty.Var(zc)));
        }
        tracksKept += slopeX.size();
        tree->Fill();
    }
 
    std::cout << TString::Format("\n%lld events read, %lld (%.1f%%) with one hit in each SiBT station\n"
                                 "target z from %.4f to %.4f cm, residuals taken at z_c = %.4f cm\n"
                                 "%lld tracks accepted (%s), %.1f per event with a beam line\n",
                                 events, eventsWithBeam, 100.0 * eventsWithBeam / std::max(events, 1LL), zStart, zEnd,
                                 zc, tracksKept, BMN_TRACKS_CUT_STR,
                                 static_cast<Double_t>(tracksKept) / std::max(eventsWithBeam, 1LL));
    if (!tracksKept)
        return -1;
 
    tree->SetEstimate(std::max(tracksKept, eventsWithBeam) + 1);
    auto quantile = [tree](const TString& expr, const TString& cut, Double_t q) {
        const Long64_t n{tree->Draw(expr, cut, "goff")};
        return n > 0 ? TMath::KOrdStat(n, tree->GetV1(), static_cast<Long64_t>(q * (n - 1))) : 0.0;
    };
    const Double_t madToSigma{1.0 / TMath::NormQuantile(0.75)};
 
    const TString vertexLike{TString::Format("(abs(slopeX - beamSlopeX) >= %g || abs(slopeY - beamSlopeY) >= %g)",
                                             BEAM_LIKE_SLOPE, BEAM_LIKE_SLOPE)};
    const Long64_t tracksVertexLike{tree->Draw("resX", vertexLike, "goff")};
    TString bothWindows{vertexLike};
 
    for (const TString p : {"X", "Y"}) {
        const TString res{"res" + p}, slope{"(slope" + p + " - beamSlope" + p + ")"}, err{"err" + p},
            beamErr{"beamErr" + p}, beamPull{"beamPull" + p};
        const Double_t chi2ndf{(p == "X" ? chi2X : chi2Y) / ndf}, chi2Other{(p == "X" ? chi2Y : chi2X) / ndf};
 
        // a fixed misalignment of the middle SiBT station shifts every pull alike,
        // so the error scale is the spread around the median, not the width around zero
        const Double_t pullOffset{quantile(beamPull, "", 0.5)},
            beamScale{madToSigma * quantile(TString::Format("abs(%s - (%g))", beamPull.Data(), pullOffset), "", 0.5)};
 
        // The target depth is a uniform spread: in the pull it enters as a variance,
        // in the window as a hard edge, so the window does not widen with WINDOW_SIGMAS
        // at large slopes. Secondaries are trimmed by shrinking the window step by step,
        // and the lines are fitted as a mean inside the window, so they follow it.
        const TString total{TString::Format("sqrt((%g * %s)^2 + %s^2)", beamScale, beamErr.Data(), err.Data())},
            depth{TString::Format("abs%s * %g", slope.Data(), halfT)}, side[2]{slope + " < 0", slope + " >= 0"};
        const Double_t slope99{quantile("abs" + slope, vertexLike, 0.99)};
        // owned by the file, which deletes them on Close()
        TProfile* const profile[2]{
            new TProfile("vertexZ" + p + "neg", "slope < 0;track slope - beam slope;beam - track at z_c, cm", 20,
                         -slope99, 0.0),
            new TProfile("vertexZ" + p + "pos", "slope > 0;track slope - beam slope;beam - track at z_c, cm", 20, 0.0,
                         slope99)};
        Double_t offset[2]{0.0, 0.0}, dz[2]{0.0, 0.0}, scale{1.0};
        Bool_t fitted[2]{kFALSE, kFALSE};
        TString window{vertexLike}, fromLine;
        for (Int_t i{0}; i < 4; i++) {
            for (Int_t s{0}; s < 2; s++) {
                tree->Draw(TString::Format("%s:%s>>%s", res.Data(), slope.Data(), profile[s]->GetName()),
                           window + " && " + side[s], "prof goff");
                // a side without tracks keeps a zero line, which selects nothing
                fitted[s] = profile[s]->GetEntries() > 0 && 0 == static_cast<Int_t>(profile[s]->Fit("pol1", "Q0"));
                if (fitted[s]) {
                    offset[s] = profile[s]->GetFunction("pol1")->GetParameter(0);
                    dz[s] = profile[s]->GetFunction("pol1")->GetParameter(1);
                }
            }
            fromLine = TString::Format("(%s - (%s) * ((%g) + (%g) * %s) - (%s) * ((%g) + (%g) * %s))", res.Data(),
                                       side[0].Data(), offset[0], dz[0], slope.Data(), side[1].Data(), offset[1], dz[1],
                                       slope.Data());
            scale = madToSigma
                    * quantile(
                        TString::Format("abs%s / sqrt(%s^2 + (%s)^2 / 3)", fromLine.Data(), total.Data(), depth.Data()),
                        window, 0.5);
            window = TString::Format("%s && abs%s < %g * %s + %s", vertexLike.Data(), fromLine.Data(),
                                     WINDOW_SIGMAS * scale, total.Data(), depth.Data());
        }
        bothWindows += " && " + window;
 
        Long64_t sideTracks[2]{}, sideInWindow[2]{};
        Double_t sideWidth[2]{};
        for (Int_t s{0}; s < 2; s++) {
            sideTracks[s] = tree->Draw(res, vertexLike + " && " + side[s], "goff");
            sideInWindow[s] = tree->Draw(res, window + " && " + side[s], "goff");
            sideWidth[s] = madToSigma * quantile("abs" + fromLine, window + " && " + side[s], 0.5);
            if (!fitted[s])
                offset[s] = dz[s] = sideWidth[s] = TMath::QuietNaN();
        }
        const Long64_t tracksInWindow{tree->Draw(res, window, "goff")};
        const Double_t totalMedian{quantile(total, window, 0.5)}, slopeMedian{quantile("abs" + slope, window, 0.5)},
            slope90{quantile("abs" + slope, window, 0.9)};
        const Long64_t thicknessWins{tree->Draw(
            res, window + TString::Format(" && (%s)^2 / 3 > (%g * %s)^2", depth.Data(), beamScale, beamErr.Data()),
            "goff")};
 
        std::cout << TString::Format("\n--- projection %s\n", p.Data());
        std::cout << TString::Format("straight line chi2/ndf over accepted tracks: %.2f\n", chi2ndf);
        if (chi2ndf > 10.0 * chi2Other)
            std::cout << "  ! far above the other projection: check the field used in reco,\n"
                         "    the hit errors and the misalignment before trusting the numbers below\n";
        std::cout << TString::Format(
            "beam at z_c: sigma %.4f cm predicted, pull scale %.2f (middle hit offset %.2f sigma)\n"
            "beam - track at z_c = offset + dz * (track slope - beam slope), fitted per sign of the slope:\n"
            "                              slope < 0   slope > 0\n"
            "  tracks                      %9lld   %9lld\n"
            "  in the vertex window, %%     %9.1f   %9.1f\n"
            "  offset, cm                  %9.4f   %9.4f\n"
            "  dz, cm                      %9.3f   %9.3f\n"
            "  width around the line, cm   %9.4f   %9.4f\n"
            "dz is the vertex z - z_c in the tracker frame: interaction depth plus the tracker z offset\n"
            "vertex window: |r - line| < %g * %.2f * total error + |track slope - beam slope| * %.4f\n"
            "total error median %.4f cm; |track slope - beam slope| in window: median %.3f, 90%% below %.3f\n"
            "thickness term (t^2/12) above beam error: %.1f%% of tracks in window\n",
            quantile(beamErr, "", 0.5), beamScale, pullOffset, sideTracks[0], sideTracks[1],
            100.0 * sideInWindow[0] / std::max(sideTracks[0], 1LL),
            100.0 * sideInWindow[1] / std::max(sideTracks[1], 1LL), offset[0], offset[1], dz[0], dz[1], sideWidth[0],
            sideWidth[1], WINDOW_SIGMAS, scale, halfT, totalMedian, slopeMedian, slope90,
            100.0 * thicknessWins / std::max(tracksInWindow, 1LL));
    }
 
    // one row per track in both windows, each carrying the count of such tracks in its event
    const Long64_t tracksWithVertex{tree->Draw("Sum$(" + bothWindows + ")", bothWindows, "goff")};
    std::cout << TString::Format(
        "\n--- both projections\n"
        "beam-like tracks (|track slope - beam slope| < %g in x and y), left out: %.1f%% of accepted\n"
        "tracks in both vertex windows: %.1f%% of accepted, %.1f per event with a beam line\n"
        "copies of the event's vertex hit, mean over the tracks that take it: %.1f\n",
        BEAM_LIKE_SLOPE, 100.0 * (tracksKept - tracksVertexLike) / tracksKept, 100.0 * tracksWithVertex / tracksKept,
        static_cast<Double_t>(tracksWithVertex) / eventsWithBeam,
        tracksWithVertex > 0 ? TMath::Mean(tracksWithVertex, tree->GetV1()) : 0.0);
 
    out.Write();
    out.Close();
    std::cout << "\n" << outPath << " : ntuple 'vertex' and the vertex z profiles\n";
    return 0;
}