review.C

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.ru
 
    Preliminary data review
*/
 
// headers below are for IDE syntax to work and they will be
// turned off during execution in ROOT macro mode
#ifndef __CLING__
#include "BmnRootIterator.h"
#include "BmnRun8cgDetModel.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TStopwatch.h"
 
#include <array>
#include <fstream>
#include <map>
#include <vector>
#endif
 
#define BMN_FILELIST "cfg/dst7466noa-full.txt"
 
#define XHIST_WIDTH 0.06
#define YHIST_WIDTH 0.6
 
// canvas and pads configuration
#define RES_X 3440
#define RES_Y 1440
#define PADS_X 4
#define PADS_Y 2
 
// filter tracks by detector element (must have hit in it)
#define DET_FILTER 201392146   // -1 to disable filtering
 
static const std::vector<Int_t> detIDs = {201588754, 201588770};
static const Double_t detMinMax[2][6] = {{-27.0, 27.0, 0.0, 16.0}, {-27.0, 27.0, -16.0, 0.0}};
 
enum
{
    hResX,
    hResY,
    hXresX,
    hYresX
};
 
std::map<Int_t, std::array<Double_t, 6>> detPos;
 
int review(int argc = 0, char** argv = nullptr)
{
    TStopwatch timer;
 
    // prepare canvas and pads
    TCanvas* pCnv = new TCanvas(   // RAM is not released in macro mode
        "canvas", "BM@N alignment");
    pCnv->SetRealAspectRatio();
    pCnv->SetFixedAspectRatio();
    pCnv->SetCanvasSize(RES_X - 20, RES_Y - 96);
    pCnv->SetWindowSize(RES_X, RES_Y);
    constexpr const Double_t padW{1.0 / PADS_X}, padH{1.0 / PADS_Y};
    TPad* pPads[PADS_Y][PADS_X];
    for (Int_t i = 0; i < PADS_Y; i++) {
        for (Int_t j = 0; j < PADS_X; j++) {
            pPads[i][j] = new TPad(Form("pad%d%d", i + 1, j + 1), Form("pad%d%d", i + 1, j + 1), padW * j, padH * i,
                                   padW * (j + 1), padH * (i + 1));
            pPads[i][j]->Draw();
        }
    }
 
    // prepare ROOT file readers and histograms
    BmnRun8cgDetModel detModel;
    detModel.Init();
    BmnRootIterator data(BMN_FILELIST, detModel);
 
    TH1D::StatOverflows(kTRUE);
    TH1* hist[PADS_Y][PADS_X];
    for (Int_t i = 0; i < PADS_Y; i++)
        for (Int_t j = 0; j < PADS_X; j++)
            hist[i][j] = nullptr;
 
    // process ROOT data
    BmnRootIterator::Hit_t detHit;
    while (!data.EndOfTracks()) {
        Bool_t filterOK{kFALSE};
        Int_t detIdx{-1};
        Double_t xSum{0.0}, ySum{0.0}, zSum{0.0}, z2Sum{0.0}, zxSum{0.0}, zySum{0.0};
        while (!data.EndOfHits()) {
            Int_t detID{data.HitDetectorID()};
            filterOK |= (detID == DET_FILTER);
            auto found{std::find(detIDs.begin(), detIDs.end(), detID)};
            if (found != detIDs.end()) {
                // hit is in one of investigated detector elements
                detHit = *data.GetHit();
                detIdx = std::distance(detIDs.begin(), found);
            } else {
                // linear regression terms
                xSum += data.HitX();
                ySum += data.HitY();
                zSum += data.HitZ();
                z2Sum += data.HitZ() * data.HitZ();
                zxSum += data.HitZ() * data.HitX();
                zySum += data.HitZ() * data.HitY();
            }
            data.NextHit();
        }
 
        if (detIdx == -1 || (!filterOK && DET_FILTER != -1)) {
            // no hits in investigated detector elements
            data.NextTrack();
            continue;
        }
        // init histograms if not yet
        if (!hist[detIdx][hResX]) {
            hist[detIdx][hResX] =
                new TH1D(Form("resX_%d", detIDs[detIdx]), "Residuals in X", 101, -XHIST_WIDTH, XHIST_WIDTH);
            hist[detIdx][hResY] =
                new TH1D(Form("resY_%d", detIDs[detIdx]), "Residuals in Y", 101, -YHIST_WIDTH, YHIST_WIDTH);
            hist[detIdx][hXresX] = new TH2D(Form("XresX_%d", detIDs[detIdx]), "Residuals in X over X", 101,
                                            detMinMax[detIdx][0], detMinMax[detIdx][1], 101, -XHIST_WIDTH, XHIST_WIDTH);
            hist[detIdx][hYresX] = new TH2D(Form("YresX_%d", detIDs[detIdx]), "Residuals in X over Y", 101,
                                            detMinMax[detIdx][2], detMinMax[detIdx][3], 101, -XHIST_WIDTH, XHIST_WIDTH);
        }
        // linear regression fit
        Double_t hitCnt{static_cast<Double_t>(data.HitsInTrack()) - 1.0}, factor{1.0 / (hitCnt * z2Sum - zSum * zSum)},
            xA{(xSum * z2Sum - zSum * zxSum) * factor}, yA{(ySum * z2Sum - zSum * zySum) * factor},
            xB{(hitCnt * zxSum - zSum * xSum) * factor}, yB{(hitCnt * zySum - zSum * ySum) * factor},
            resX{detHit.GetX() - xA - xB * detHit.GetZ()}, resY{detHit.GetY() - yA - yB * detHit.GetZ()};
        // residuals in target detector
        hist[detIdx][hResX]->Fill(resX);
        hist[detIdx][hResY]->Fill(resY);
        ((TH2D*)hist[detIdx][hXresX])->Fill(detHit.GetX(), resX);
        ((TH2D*)hist[detIdx][hYresX])->Fill(detHit.GetY(), resX);
 
        data.NextTrack();
    }
 
    std::ofstream f("positions.csv", std::ofstream::out);
    for (auto i = detPos.begin(); i != detPos.end(); i++) {
        std::cout << i->first << ",";
        f << i->first << ",";
        for (Int_t j = 0; j < 6; j++) {
            std::cout << i->second[j] << ",";
            f << i->second[j] << ",";
        }
        std::cout << "\n";
        f << "\n";
    }
    f.close();
 
    // draw histograms
    TH1D* h1;
    TH2D* h2;
    for (Int_t i = 0; i < PADS_Y; i++) {
        if (!hist[i][0])
            continue;
        pPads[i][hResX]->cd();
        h1 = (TH1D*)hist[i][hResX];
        h1->SetFillColor(kCyan);
        h1->Draw("HIST");
        pPads[i][hResY]->cd();
        h1 = (TH1D*)hist[i][hResY];
        h1->SetFillColor(kTeal);
        h1->Draw("HIST");
        pPads[i][hXresX]->cd();
        h2 = (TH2D*)hist[i][hXresX];
        h2->SetStats(1);
        h2->Draw("CONTZ");
        pPads[i][hYresX]->cd();
        h2 = (TH2D*)hist[i][hYresX];
        h2->SetStats(1);
        h2->Draw("CONTZ");
    }
 
    // no memory cleanup in ROOT macro mode
    std::cout << "Done in " << timer.RealTime() << "s" << std::endl;
    return 0;
}