compareA.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
 
    Aligmnment quality evaluation macro
*/
 
// headers below are for IDE syntax to work and they will be
// turned off during execution in ROOT macro mode
#ifndef __CLING__
#include "BmnDataReader.h"
#include "CbmStsTrack.h"
#include "CbmVertex.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TPad.h"
#include "TPaveText.h"
#include "TString.h"
#include "TTreeReaderArray.h"
#include "TTreeReaderValue.h"
 
#include <iomanip>
#include <iostream>
#endif
 
#define BMN_FILELIST_NA "cfg/dst_mini_fe.txt"
#define BMN_FILELIST_A "cfg/dst_mini_re.txt"
constexpr static const char* BMN_RESULT_PDF{
    nullptr};   //{"/home/igor/DATA/bmn/dst/25.04.0/compare/cmp_ali-noa_7466_ev0_p4.pdf"};
 
// canvas and pads configuration
#define RES_X 3440
#define RES_Y 1440
#define PADS_X 4
#define PADS_Y 2
 
enum eAligned
{   // aligned or not
    kNA = 0,
    kAL = 1
};
constexpr Int_t kDEF{kAL}, kSEC{kNA};
 
enum eCmpHists
{   // comparison histograms enumeration
    kTracksCount,
    kVertexCount,
    kTracksChi2perNDF,
    kVertexX,
    kVertexY,
    kVertexZ,
    kVertexXYbefore,
    kVertexXYafter,
    kCount
};
 
// historgam positioning on pads
constexpr const Int_t histPos[eCmpHists::kCount][2]{
    {1, 1},   // kTracksCount
    {1, 0},   // kVertexCount
    {1, 2},   // kTracksChi2perNDF
    {0, 0},   // kVertexX
    {0, 1},   // kVertexY
    {0, 2},   // kVertexZ
    {1, 3},   // kVertexXYbefore
    {0, 3}    // kVertexXYafter
};
 
int compareA(int argc = 0, char** argv = nullptr)
{
    // 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 (UInt_t i = 0; i < PADS_Y; i++) {
        for (UInt_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
    TH1* pHist[2][eCmpHists::kCount];
    BmnDataReader* pData[2];
    TTreeReaderArray<CbmStsTrack>* pTracks[2];
    TTreeReaderValue<CbmVertex>* pVertex[2];
    for (Int_t al = kNA; al <= kAL; al++) {
        pData[al] = new BmnDataReader(al ? BMN_FILELIST_A : BMN_FILELIST_NA);
        pTracks[al] = new TTreeReaderArray<CbmStsTrack>(pData[al]->Reader(), "StsVector");
        pVertex[al] = new TTreeReaderValue<CbmVertex>(pData[al]->Reader(), "MpdVertex.");
        pHist[al][kVertexCount] = new TH1D(al ? "VC-AL" : "VC-NA", "Vertex count (tracks > 1)", 4, -1.0, 2.0);
        pHist[al][kTracksChi2perNDF] =
            new TH1D(al ? "TChi-AL" : "TChi-NA", "Tracks Chi2 per NDF (cut = 5.0)", 100, 0.0, 5.0);
        pHist[al][kVertexX] = new TH1D(al ? "VX-AL" : "VX-NA", "Vertex X distribution", 100, -6.0, 6.0);
        pHist[al][kVertexY] = new TH1D(al ? "VY-AL" : "VY-NA", "Vertex Y distribution", 100, -6.0, 6.0);
        pHist[al][kVertexZ] = new TH1D(al ? "VZ-AL" : "VZ-NA", "Vertex Z distribution", 100, -2.0, 2.0);
    }
    pHist[kDEF][kVertexXYbefore] = new TH2D("VXY-NA", "Non-aligned vertex XY", 100, -2.0, 2.0, 100, -2.0, 2.0);
    pHist[kDEF][kVertexXYafter] = new TH2D("VXY-AL", "Aligned vertex XY", 100, -2.0, 2.0, 100, -2.0, 2.0);
    pHist[kDEF][kTracksCount] = new TH1D("TC", "Tracks count (hits > 1)", 4, -1.0, 2.0);
    pHist[kSEC][kVertexXYbefore] = nullptr;
    pHist[kSEC][kVertexXYafter] = nullptr;
    pHist[kSEC][kTracksCount] = nullptr;
 
    // fill histograms
    TH1D::StatOverflows(kFALSE);   // need to ignore overflows, otherwise
                                   // vertex Z histogram may show increase of StdDev
    Int_t vertexTotal[2]{0, 0}, tracksTotal[2]{0, 0};
    for (Int_t al = kNA; al <= kAL; al++) {
        Double_t dCounter{static_cast<Double_t>(al)};
        TH2D* pHistVertexXY = (TH2D*)(al == kAL ? pHist[kDEF][kVertexXYafter] : pHist[kDEF][kVertexXYbefore]);
        while (pData[al]->Reader().Next()) {
            // vertex data
            CbmVertex& vertex = **pVertex[al];
            if (vertex.GetNTracks() < 2)   // ignore vertices with less than 2 tracks
                continue;
            vertexTotal[al]++;
            pHist[kDEF][kVertexCount]->Fill(dCounter);
            pHist[al][kVertexX]->Fill(vertex.GetX());
            pHist[al][kVertexY]->Fill(vertex.GetY());
            pHist[al][kVertexZ]->Fill(vertex.GetZ());
            pHistVertexXY->Fill(vertex.GetX(), vertex.GetY());
            // tracks data
            for (CbmStsTrack& track : *pTracks[al]) {
                if (track.GetNStsHits() < 2)   // ignore tracks with less than 2 hits
                    continue;
                tracksTotal[al]++;
                pHist[kDEF][kTracksCount]->Fill(dCounter);
                pHist[al][kTracksChi2perNDF]->Fill(track.GetChi2() / track.GetNDF());
            }
        }
    }
 
    // rescale chi2/hit histogram
    Double_t tracksIncrease{static_cast<Double_t>(tracksTotal[kAL]) / tracksTotal[kNA]};
    pHist[kNA][kTracksChi2perNDF]->Scale(tracksIncrease);
 
    // summarize results
    Double_t rTracksIncrease{100.0 * tracksIncrease - 100.0},
        rVertexIncrease{100.0 * vertexTotal[kAL] / vertexTotal[kNA] - 100.0},
        rTracksChi2perNDFMeanDecrease{
            100.0 * pHist[kAL][kTracksChi2perNDF]->GetMean() / pHist[kNA][kTracksChi2perNDF]->GetMean() - 100.0},
        rVertexXStdDevDecrease{100.0 * pHist[kAL][kVertexX]->GetStdDev() / pHist[kNA][kVertexY]->GetStdDev() - 100.0},
        rVertexYStdDevDecrease{100.0 * pHist[kAL][kVertexY]->GetStdDev() / pHist[kNA][kVertexY]->GetStdDev() - 100.0},
        rVertexZStdDevDecrease{100.0 * pHist[kAL][kVertexZ]->GetStdDev() / pHist[kNA][kVertexZ]->GetStdDev() - 100.0};
 
    // report to stdout
    using namespace std;
    cout << fixed << setprecision(1);
    cout << "Tracks count increase ......................... " << rTracksIncrease << "%\n";
    cout << "Vertex count increase ......................... " << rVertexIncrease << "%\n";
    cout << "Mean tracks Chi2 per NDF decrease ............. " << rTracksChi2perNDFMeanDecrease << "%\n";
    cout << "Vertex X distribution StdDev decrease ......... " << rVertexXStdDevDecrease << "%\n";
    cout << "Vertex Y distribution StdDev decrease ......... " << rVertexYStdDevDecrease << "%\n";
    cout << "Vertex Z distribution StdDev decrease ......... " << rVertexZStdDevDecrease << "%\n";
    cout << "Vertex Z distribution StdDev before / after ..... " << setprecision(3) << pHist[kNA][kVertexZ]->GetStdDev()
         << " / " << pHist[kAL][kVertexZ]->GetStdDev() << "\n";
 
    // report to GUI
    TH1* hst[2];
    TPaveText* pave;
    TString hstOpt[2];
    EColor colorNA;
    Double_t value;
    for (Int_t hstEN = 0; hstEN < eCmpHists::kCount; hstEN++) {
        pPads[histPos[hstEN][0]][histPos[hstEN][1]]->cd();
        hst[kNA] = pHist[kNA][hstEN];
        hst[kAL] = pHist[kAL][hstEN];
        // default histogram configuration
        hst[kDEF]->SetLineColor(kRed);
        hst[kDEF]->SetLineWidth(2);
        hst[kDEF]->SetStats(0);
        if (hst[kSEC])
            hst[kSEC]->SetStats(0);
        hstOpt[kNA] = "HIST";
        hstOpt[kAL] = "C SAME";
        // every TPaveText object is left in heap due to ROOT macro drawing specifics
        pave = new TPaveText(0.6, 0.85, 0.95, 0.92, "BR NDC");
        switch (hstEN) {
            case kTracksCount:
            case kVertexCount:
                hst[kDEF]->GetXaxis()->SetBinLabel(1, "");
                hst[kDEF]->GetXaxis()->SetBinLabel(kNA + 2, "Not aligned");
                hst[kDEF]->GetXaxis()->SetBinLabel(kAL + 2, "Aligned");
                hst[kDEF]->GetXaxis()->SetBinLabel(4, "");
                hst[kDEF]->SetFillColor(hstEN == kTracksCount ? kCyan : kTeal);
                hst[kDEF]->SetBarWidth(0.8);
                hst[kDEF]->SetBarOffset(0.1);
                hstOpt[kDEF] = "BAR";
                pave->AddText(hstEN == kTracksCount ? Form("Tracks increase: %+5.1f%%", rTracksIncrease)
                                                    : Form("Vertex increase: %+5.1f%%", rVertexIncrease));
                break;
            case kTracksChi2perNDF:
                colorNA = kCyan;
                pave->AddText(Form("Mean decrease: %+5.1f%%", rTracksChi2perNDFMeanDecrease));
                break;
            case kVertexX:
            case kVertexY:
            case kVertexZ:
                colorNA = (hstEN == kVertexZ) ? kGreen : kTeal;
                if (hstEN == kVertexX)
                    value = rVertexXStdDevDecrease;
                else if (hstEN == kVertexY)
                    value = rVertexYStdDevDecrease;
                else if (hstEN == kVertexZ)
                    value = rVertexZStdDevDecrease;
                pave->AddText(Form("StdDev decrease: %+5.1f%%", value));
                break;
            case kVertexXYbefore:
            case kVertexXYafter:
                hstOpt[kDEF] = "CONTZ";
                delete pave, pave = nullptr;
                break;
            default:
                hst[kNA] = nullptr;
                hst[kAL] = nullptr;
                delete pave, pave = nullptr;
        }
        if (hst[kNA]) {
            hst[kNA]->SetFillColor(colorNA);
            hst[kNA]->SetLineColor(colorNA);
            if (hstOpt[kNA] != "CONTZ")
                hst[kNA]->GetYaxis()->SetRangeUser(0.0, hst[kAL]->GetMaximum() * 1.1);
            hst[kNA]->Draw(hstOpt[kNA]);
        }
        if (hst[kAL])
            hst[kAL]->Draw(hstOpt[kAL]);
        if (pave)
            pave->Draw();
    }
 
    // report to PDF
    if (BMN_RESULT_PDF)
        pCnv->SaveAs(BMN_RESULT_PDF);
 
    // partial cleanup
    for (UInt_t i = 0; i < 2; i++) {
        delete pData[i];
        delete pTracks[i];
        delete pVertex[i];
    }
 
    return 0;
}