BmnAlignerApply.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.ru
 
    BmnAlignerApply class implementation
*/
 
#include "BmnAlignerApply.h"
 
#include "BmnGemStripStationSet.h"
#include "BmnSiliconStationSet.h"
#include "DstEventHeader.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "TObjString.h"
#include "TSystem.h"
#include "fairlogger/Logger.h"
 
#include <iostream>
 
InitStatus BmnAlignerApply::Init()
{
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        LOG(error) << " <!> BmnAligner::Init() : FairRootManager not found";
        return kERROR;
    }
    fpStsHitsArray = (TClonesArray*)ioman->GetObject("StsHit");
    if (!fpStsHitsArray) {
        LOG(error) << " <!> BmnAligner::Init() : Failed to get StsHit array";
        return kERROR;
    }
    try {
        if (!LoadSolution(fAlignFilePath, kFALSE, kFALSE)) {
            LOG(error) << " <!> BmnAlignerApply::Init() : Failed to load alignment solution from file "
                       << fAlignFilePath;
            return kERROR;
        }
    } catch (std::exception& e) {
        LOG(error) << " <!> BmnAlignerApply::Init() : " << e.what();
        return kERROR;
    } catch (...) {
        LOG(error) << " <!> BmnAlignerApply::Init() : Unknown exception";
        return kERROR;
    }
    if (!fLorentzCorrectionApply)
        return kSUCCESS;
 
    /*
     * Lorentz corrections section-----------------------------------------------------------------
     *
     * The code below is directly taken from BmnToCbmHitConverter::ApplyAlignment(CbmStsHit*).
     * It implements the Lorentz correction of hit positions and was not revised, except few
     * lines marked with tag PLV. Indentation is the kept as in the original code.
     */
 
    FairField* mField = FairRunAna::Instance()->GetField();
    fieldFlag = (TMath::Abs(mField->GetBy(0, 0, 135)) > 0.5) ? kTRUE : kFALSE;
    TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
    // PLV : begin
    // fGemConfigFile and fSilConfigFile initialization moved from ctor
    // original code:
    // , fGemConfigFile("GemRun8.xml")
    // , fSilConfigFile("SiliconRun8.xml")
    // modified code:
    TString fGemConfigFile("GemRun8.xml"), fSilConfigFile("SiliconRun8.xml");
    // BmnGemStripStationSet and BmnSiliconStationSet moved from class fields to local variables
    // original code:
    // GemStationSet = new BmnGemStripStationSet(gPathConfig + "/parameters/gem/XMLConfigs/" + fGemConfigFile);
    // SilStationSet = new BmnSiliconStationSet(gPathConfig + "/parameters/silicon/XMLConfigs/" + fSilConfigFile);
    // modified code:
    BmnGemStripStationSet* GemStationSet =
        new BmnGemStripStationSet(gPathConfig + "/parameters/gem/XMLConfigs/" + fGemConfigFile);
    BmnSiliconStationSet* SilStationSet =
        new BmnSiliconStationSet(gPathConfig + "/parameters/silicon/XMLConfigs/" + fSilConfigFile);
    // PLV : end
 
    static int irun = -1;
    // PLV : begin
    // zStat is not used in calculations
    // static Double_t zStat[19] = {0};
    // PLV : end
 
    if (irun <= 0)
        irun = ((DstEventHeader*)FairRootManager::Instance()->GetObject("DstEventHeader."))->GetRunId();   // AZ-010723
 
    // Lorentz shifts
    // AZ-010723 if (hLorCor[0][0] == nullptr) {
    if (irun <= 0)
        // PLV : begin
        // all the section moved from method ApplyAlignment() - return type is different
        // original code:
        // return;
        // modified code:
        return kERROR;
    // PLV : end
    else if (hLorCorX[0][0] == nullptr) {
        // Station Z-positions
        int nsta = SilStationSet->GetNStations();
        // PLV : begin
        // zStat is not used in calculations
        // for (int ista = 0; ista < nsta; ++ista)
        //     zStat[ista] = SilStationSet->GetZStationPosition(ista);
        // PLV : end
        int nsta1 = GemStationSet->GetNStations();
        // PLV : begin
        // zStat is not used in calculations
        // for (int ista = 0; ista < nsta1; ++ista)
        //     zStat[ista + nsta] = GemStationSet->GetZStationPosition(ista);
        // PLV : end
        // PLV : begin
        // non-necessary console output removed
        // if (fVerbose > 0)
        //     for (int j = 0; j < nsta + nsta1; ++j)
        //         cout << " Stat: " << j << " " << zStat[j] << endl;
        // PLV : end
 
        map<int, vector<Double_t>> lorCorX[2], lorCorY[2];
        int nx, ny;
 
        for (int igem = 0; igem < 2; ++igem) {
            ReadCorrections(irun, igem, lorCorX, lorCorY, nx, ny);
            Double_t dx = 160.0 / nx, dy = 45.0 / ny;
            // Double_t dx = 20.0, dy = 15.0;
 
            for (int iupd = 0; iupd < 2; ++iupd) {
                int istb = (igem == 0) ? 0 : nsta;
                int iste = (igem == 0) ? nsta : nsta + nsta1;
                if (iste > (int)lorCorX[0].size())
                    --iste;   // exclude last small GEM
 
                for (int ist = istb; ist < iste; ++ist) {
                    TString hName = "hLorCorX_";
                    hName += iupd;
                    hName += "_";
                    hName += ist;
                    // hLorCor[iupd][ist] = new TH2D(hName, "", 8, -80, 80, 3, 0, 45);
                    int nxfsd = 0;
                    if (igem)
                        hLorCorX[iupd][ist] = new TH2D(hName, "", nx, -80, 80, ny, 0, 45);
                    else {
                        nxfsd = lorCorX[0][ist + 1].size();
                        hLorCorX[iupd][ist] = new TH2D(hName, "", nxfsd, 0, nxfsd, ny, 0, ny);
                    }
 
                    hName = "hLorCorY_";
                    hName += iupd;
                    hName += "_";
                    hName += ist;
                    if (igem)
                        hLorCorY[iupd][ist] = new TH2D(hName, "", nx, -80, 80, ny, 0, 45);
                    else
                        hLorCorY[iupd][ist] = new TH2D(hName, "", nxfsd, 0, nxfsd, ny, 0, ny);
 
                    if (igem) {
                        // for (int j = 0; j < nx*ny; ++j) {
                        for (int j = 0; j < 24; ++j) {
                            // int ireg1 = j - 13;
                            // if (ireg1 < 0) ++ireg1;
                            int ireg1 = j - 12;
                            int ireg = TMath::Abs(ireg1);
                            if (ireg1 < 0)
                                --ireg;
                            int iy2 = ireg / 4;
                            int ix2 = ireg % 4;
                            Double_t x = (ix2 + 0.5) * dx;
                            if (ireg1 < 0)
                                x = -x;
                            Double_t y = (iy2 + 0.5) * dy;
                            hLorCorX[iupd][ist]->Fill(x, y, lorCorX[iupd][ist + 1][j]);   // X-corrections
                            if (lorCorY[iupd].find(ist + 1) != lorCorY[iupd].end())
                                hLorCorY[iupd][ist]->Fill(x, y, lorCorY[iupd][ist + 1][j]);   // Y-corrections
                        }   // for (int j = 0;
                    } else {
                        for (int j = 0; j < nxfsd; ++j) {
                            hLorCorX[iupd][ist]->Fill(j + 0.5, 0.5, lorCorX[iupd][ist + 1][j]);   // X-corrections
                            hLorCorY[iupd][ist]->Fill(j + 0.5, 0.5, lorCorY[iupd][ist + 1][j]);   // Y-corrections
                        }
                    }   // if (igem)
                }   // for (int ist = istb;
            }   // for (int iupd = 0;
        }   // for (int igem = 0;
    }   // else if (hLorCorX[0][0] == nullptr)
 
    /*
     * End of Lorentz corrections section ---------------------------------------------------------
     */
    return kSUCCESS;
}
 
void BmnAlignerApply::Exec(Option_t* opt)
{
    for (Int_t i = 0, n = fpStsHitsArray->GetEntriesFast(); i < n; i++) {
        CbmStsHit* hit = static_cast<CbmStsHit*>(fpStsHitsArray->At(i));
        SVectGL& align = fpResultCurrent->A(fpDetModel->IdxFromHit(hit->GetDetectorID()));
        // keep alignment application code inside BmnMeasureModel
        // to prevent inconsistency in sign selection and rotation sequence
        fpMeasureModel->AlignHit(*hit, align);
        if (!fLorentzCorrectionApply)
            continue;
 
        /*
         * Lorentz corrections section-----------------------------------------------------------------
         *
         * The code below is directly taken from BmnToCbmHitConverter::ApplyAlignment(CbmStsHit*).
         * It implements the Lorentz correction of hit positions and was not revised, except few
         * lines marked with tag PLV. Indentation is the kept as in the original code.
         */
 
        int ista = hit->GetStationNr() - 1;
        int isec = hit->GetSectorNr() - 1;
        if (fieldFlag && hit->GetSystemId() == kGEM) {
            // Lorentz corrections
            int iupd = isec / 4;
            Double_t xh = hit->GetX(), yh = TMath::Abs(hit->GetY()), xh1 = xh, yh1 = yh;
            TH2D* h2 = hLorCorX[iupd][ista];
            if (h2 == nullptr)
                // PLV : begin
                // in case of non-initialized h2 we must proceed to the next hit
                // original code:
                // return;
                // modified code:
                continue;
            // PLV : end
            if (xh1 >= h2->GetXaxis()->GetXmax())
                xh1 = h2->GetXaxis()->GetXmax() - 0.1;
            else if (xh1 < h2->GetXaxis()->GetXmin())
                xh1 = h2->GetXaxis()->GetXmin() + 0.1;
            if (yh1 >= h2->GetYaxis()->GetXmax())
                yh1 = h2->GetYaxis()->GetXmax() - 0.1;
            else if (yh1 < h2->GetYaxis()->GetXmin())
                yh1 = h2->GetYaxis()->GetXmin() + 0.1;
            // if (yh < h2->GetYaxis()->GetXmin() || yh >= h2->GetYaxis()->GetXmax()) cout << " yy " << isec << " " <<
            // yh << endl;
            Double_t dx = h2->Interpolate(xh1, yh1);
            // dx = 0; //AZ-060723 - no corrections
            hit->SetX(xh - dx * 0.65);
            h2 = hLorCorY[iupd][ista];
            Double_t dy = h2->Interpolate(xh1, yh1);
            if (iupd)
                dy = -dy;
            // cout << iupd << " " << ista << " " << dy << endl;
            Double_t ynew = yh - dy * 0.65;
            if (iupd)
                ynew = -ynew;
            hit->SetY(ynew);
        } else if (hit->GetSystemId() != kGEM) {
            int iupd = isec / hLorCorX[0][ista]->GetNbinsX();
            int ireg = isec % hLorCorX[0][ista]->GetNbinsX();
            TH2D* h2 = hLorCorX[iupd][ista];
            Double_t dx = h2->GetBinContent(ireg + 1, 1);
            hit->SetX(hit->GetX() - dx * 0.65);
            h2 = hLorCorY[iupd][ista];
            Double_t dy = h2->GetBinContent(ireg + 1, 1);
            hit->SetY(hit->GetY() - dy * 0.65);
        }
 
        /*
         * End of Lorentz corrections section ---------------------------------------------------------
         */
 
    }   // end of loop over hits
}
 
void BmnAlignerApply::Finish()
{
    // we could do hLorCorX and hLorCorY cleanup here but
    // it would interfere with ROOT logic of working with histograms
}
 
/*
 * Lorentz corrections section-----------------------------------------------------------------
 *
 * The code below is directly taken from BmnToCbmHitConverter::ApplyAlignment(CbmStsHit*).
 * It implements the Lorentz correction of hit positions and was not revised, except few
 * lines marked with tag PLV. Indentation is the kept as in the original code.
 */
 
// PLV : begin
// method ReadCorrections() header modified for BmnAlignerApply class
// original code:
// void BmnToCbmHitConverter::ReadCorrections(int irun,
//                                            int igem,
//                                            std::map<int, std::vector<Double_t>>* lorCorX,
//                                            std::map<int, std::vector<Double_t>>* lorCorY,
//                                            int& nx,
//                                            int& ny)
// modified code:
void BmnAlignerApply::ReadCorrections(int irun,
                                      int igem,
                                      std::map<int, std::vector<Double_t>>* lorCorX,
                                      std::map<int, std::vector<Double_t>>* lorCorY,
                                      int& nx,
                                      int& ny)
// PLV : end
{
    // Read file with "effective" Lorentz-corrections
 
    ifstream fin;
    TString fileName, chline;
 
    fileName = gSystem->ExpandPathName("$VMCWORKDIR");
    // AZ-030825 fileName += "/input/";
    fileName += "/parameters/reco/";   // AZ-030825
    if (igem)
        fileName += "LorentzGemCorr.txt";
    else
        fileName += "LorentzFsdCorr.txt";
    fin.open(fileName.Data());
    if (!fin.is_open()) {   // AZ-040825
        cout << " !!! Lorentz correction file " << fileName << " not found !!!" << endl;
        exit(1);
    }
 
    map<int, std::size_t> runsMap;
    TObjArray* tokens = nullptr;
 
    while (1) {
        // while (!fin.eof()) {
        size_t pos = fin.tellg();
        chline.ReadLine(fin);
        if (fin.eof())
            break;
        if (chline.Contains("Run")) {
            // new run
            // cout << chline << endl;
            tokens = chline.Tokenize(" ");
            TString run = ((TObjString*)tokens->UncheckedAt(1))->String();
            run.Atoi();
            // runsMap[run.Atoi()] = fin.tellg();
            runsMap[run.Atoi()] = pos;
            // cout << run.Atoi() << " " << pos << endl;
        }
        if (tokens) {
            tokens->Delete();
            delete tokens;
            tokens = nullptr;
        }
    }
    // if (chline[0] == '*') continue; // comment
    // cout << runsMap.size() << endl;
 
    auto iter = runsMap.upper_bound(irun);
    if (iter != runsMap.begin())
        --iter;
    // PLV : begin
    // non-necessary console output removed
    // cout << " ---> Using corrections for run " << iter->first << " " << iter->second << " " << igem << endl;
    // PLV : end
 
    fin.clear();
    fin.seekg(iter->second);
    chline.ReadLine(fin);
    // cout << chline << endl;
 
    // Read detectors
    while (1) {
        chline.ReadLine(fin);
        if (fin.eof())
            break;
        if (chline.Contains("Station")) {
            // cout << chline << endl;
            tokens = chline.Tokenize(" ");
            TString stat = ((TObjString*)tokens->UncheckedAt(1))->String();
            int ista = stat.Atoi();
            int iupd = 0;   // half-station
            if (chline.Contains("lower"))
                iupd = 1;
            nx = ((TObjString*)tokens->UncheckedAt(3))->String().Atoi();
            ny = ((TObjString*)tokens->UncheckedAt(4))->String().Atoi();
            int ixdata = 0;
            if (chline.Contains("X-"))
                ixdata = 1;
            // cout << ista << " " << iupd << " " << nx << " " << ny << endl;
            if (nx * ny < 1)
                cout << " !!! Wrong input line !!! " << chline << " " << nx << " " << ny << endl;
            if (tokens) {
                tokens->Delete();
                delete tokens;
                tokens = nullptr;
            }
            chline.ReadLine(fin);
            tokens = chline.Tokenize(", ");
            int nreg = tokens->GetEntriesFast();
            if (nreg != nx * ny) {
                cout << " !!! Wrong input line !!! " << chline << " " << nreg << " " << nx * ny << endl;
                exit(0);
            }
            auto lorCor = vector<Double_t>(nreg);
            // cout << tokens->GetEntriesFast() << " " << ((TObjString*)tokens->UncheckedAt(0))->String() << endl;
            for (int itok = 0; itok < nreg; ++itok) {
                Double_t corr = ((TObjString*)tokens->UncheckedAt(itok))->String().Atof();
                // lorCor[iupd][ista][itok] = corr;
                lorCor[itok] = corr;
            }
            if (ixdata)
                lorCorX[iupd][ista] = lorCor;
            else
                lorCorY[iupd][ista] = lorCor;
            if (tokens) {
                tokens->Delete();
                delete tokens;
                tokens = nullptr;
            }
        } else if (chline.Contains("*")) {
            continue;   // comment
        } else
            break;
    }   // while (1)
    fin.close();
}
 
/*
 * End of Lorentz corrections section ---------------------------------------------------------
 */