BmnStsVectorFinderV9.cxx

Go to the documentation of this file.

// -------------------------------------------------------------------------
// -----                  BmnStsVectorFinder source file               -----
// -----  Created 3/12/21  by A.Zinchenko, D.Zinchenko, R.Zinchenko    -----
// -------------------------------------------------------------------------
#include "BmnStsVectorFinderV9.h"
 
#include "BmnStsMatchTracks.h"
#include "CbmKF.h"
#include "CbmKFTrack.h"
#include "CbmStsCluster.h"
#include "CbmStsDigiMatch.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsHit.h"
#include "CbmStsPoint.h"
#include "CbmStsSector.h"
#include "CbmStsStation.h"
#include "CbmStsTrack.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "TMVA/Tools.h"
 
#include <TClonesArray.h>
#include <TF1.h>
#include <TGraph.h>
#include <TMath.h>
#include <TObjString.h>
#include <iomanip>
#include <iostream>
#include <map>
#include <set>
#include <vector>
 
using std::cout;
using std::endl;
using std::map;
using std::multimap;
using std::pair;
using std::set;
using std::vector;
 
#include <TStopwatch.h>
 
// std::ofstream outf("log.txt");
 
static Double_t workTime = 0.0;
 
// -----   Default constructor   -------------------------------------------
BmnStsVectorFinderV9::BmnStsVectorFinderV9()
    : FairTask("STS Vector Finder V9")
    , fHitArray(NULL)
    , fTrackArray(NULL)
    , fEvent(0)
    , fPass(0)
    , fNsta(0)
    ,
    // fNbranches(5),   // AZ-021125 - for tracking in windows
    fNbranches(7)
    ,   // AZ-230625
    // fNbranches(10), // AZ-270625
    fExact(0)
    , fExactSel(-9)
    , discarded(0)
    , fOut(0)   // 0,2,3,6
                // AZ-240626 fNhitsMin(nullptr)
    , fMatBudgetFileName("")
    , fUseTMVA(0)
    ,   // 0,2,3 //AZ-300922
    fTMVAtracks(1)
    ,   // AZ-040525 - discriminate good from fake tracks
    fNSi(4)
    ,   // AZ-051022
    fGemConfigFile("GemRun8.xml")
    , fSilConfigFile("SiliconRun8.xml")
    , fOutBranch("StsVector")
    , fKillHitsTxP(0)
    , fTracksToKill(nullptr)
    , fHitConverter(nullptr)
{
    for (int j = 0; j < 4; ++j)
        fClusArray[j] = nullptr;
    // fOutf = new ofstream("log.txt"); // AZ-290525 - for auxiliary output
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
BmnStsVectorFinderV9::~BmnStsVectorFinderV9() {}
// -------------------------------------------------------------------------
 
// -----   Public method Init   --------------------------------------------
InitStatus BmnStsVectorFinderV9::Init()
{
 
    // Get RootManager
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        cout << "-E- BmnStsVectorFinderV9::Init: " << "RootManager not instantiated!" << endl;
        return kFATAL;
    }
 
    // Get input arrays
    fClusArray[0] = (TClonesArray*)ioman->GetObject("StsCluster");
    if (!fClusArray[0]) {
        // cout << "-W- BmnStsVectorFinderV9::Init: "
        //   << "No StsCluster array!" << endl;
        // return kERROR;
        fClusArray[0] = (TClonesArray*)ioman->GetObject("BmnSiliconUpperCluster");
        fClusArray[1] = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
        fClusArray[2] = (TClonesArray*)ioman->GetObject("BmnGemUpperCluster");
        fClusArray[3] = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
        for (int j = 0; j < 4; ++j) {
            // AZ-241222 if ( ! fClusArray[j] ) {
            if (!fClusArray[0] && !fClusArray[2]) {   // AZ-241222
                cout << "-W- BmnStsVectorFinderV9::Init: " << "No BmnCluster array!" << j << endl;
                return kERROR;
            }
        }
    }
 
    fHitArray = (TClonesArray*)ioman->GetObject("StsHit");
    if (!fHitArray) {
        // cout << "-W- BmnStsVectorFinderV9::Init: "
        //   << "No StsHit array!" << endl;
        return kERROR;
    }
 
    fTrackArray = (TClonesArray*)ioman->GetObject("StsTrack");
    if (!fTrackArray) {
        // cout << "-W- BmnStsVectorFinderV9::Init: "
        //   << "No StsTrack array!" << endl;
        // AZ-100722 return kERROR;
    }
 
    fDigiMatches = (TClonesArray*)ioman->GetObject("StsDigiMatch");
    if (!fDigiMatches) {
        // cout << "-W- BmnStsVectorFinderV9::Init: "
        //   << "No StsDigiMatch array!" << endl;
        // return kERROR;
    }
 
    fStsPoints = (TClonesArray*)ioman->GetObject("StsPoint");
    if (!fStsPoints) {
        // cout << "-W- BmnStsVectorFinderV9::Init: "
        //   << "No StsPoint array!" << endl;
        // AZ-181222 return kERROR;
    }
 
    fSilPoints = (TClonesArray*)ioman->GetObject("SiliconPoint");
    if (!fSilPoints) {
        // cout << "-E- BmnStsVectorFinderV9::Init: No SiliconPoint array!" << endl;
    }
 
    fVspPoints = (TClonesArray*)ioman->GetObject("VSPPoint");
    if (!fVspPoints) {
        // cout << "-E- BmnStsVectorFinder::Init: No SiliconPoint array!" << endl;
    }
 
    // Create and register output array
    fVectorArray = new TClonesArray("CbmStsTrack");
    // AZ-041125 ioman->Register("StsVector", "STS", fVectorArray, kTRUE);
    ioman->Register(fOutBranch, "STS", fVectorArray, kTRUE);   // AZ-041125
 
    if (fOut / 3 > 0) {   // AZ-260922
        // Store triplets
        // AZ-240623 fMap3OutPtr = &fMap3Out;
        // AZ-240623 ioman->RegisterAny("Triplets", fMap3OutPtr, kTRUE);
        fStream3out.open("triplets.txt", std::ofstream::out);
    }
    if (fOut && fOut % 2 == 0) {   // AZ-260922
        // Store doublets
        // AZ-240623 fMap2OutPtr = &fMap2Out;
        // AZ-240623 ioman->RegisterAny("Doublets", fMap2OutPtr, kTRUE);
        fStream2out.open("doublets.txt", std::ofstream::out);
    }
 
    fXyzv[0] = 0.4;
    fXyzv[1] = 0.1;
    fXyzv[2] = 0.0;
 
    fitter.Init();
 
    // Define logic
    vector<int> nHitsMin = {7, 5, 4, 4};   // min number of hits per track vs iteration - Run 7
    // static Int_t nHitsMin[20] =   {7,    5,    4,    3}; // min number of hits per track vs iteration - Run 7
    // static Double_t dTanX[20] =   {0.001, 0.003, 0.002, 0.008}; // window size in TanX vs iteration
    vector<Double_t> dTanX = {0.002, 0.003, 0.002, 0.008};        // window size in TanX vs iteration
    vector<Double_t> dTanXTMVA = {0.003, 0.0045, 0.003, 0.012};   // window size in TanX vs iteration - for TMVA
    vector<Double_t> dTanXB0 = {0.001, 0.002, 0.002, 0.002};      // window size in TanX vs iteration - for zero field
    // static Double_t dTanY[20] =   {0.01, 0.02, 0.01, 0.02}; // window size in TanY vs iteration
    // AZ-230623 static Double_t dTanY[20] =   {0.01, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
    vector<Double_t> dTanY = {0.02, 0.02, 0.02, 0.03};   // window size in TanY vs iteration
    // static Double_t dTanY[20] = {0.05, 0.05, 0.05, 0.1};         // window size in TanY vs iteration - 221225 - test
    vector<Double_t> dTanYTMVA = {0.015, 0.03, 0.03, 0.045};   // window size in TanY vs iteration - for TMVA
    vector<Double_t> dTanYB0 = {0.01, 0.02, 0.02, 0.02};       // window size in TanY vs iteration - for zero field
    // AZ-220623 static Double_t dTanY3[20] =  {0.01, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
    vector<Double_t> dTanY3 = {0.03, 0.02, 0.015, 0.015};       // AZ-220623 window size in TanY vs iteration
    vector<Double_t> dTanY3TMVA = {0.015, 0.03, 0.03, 0.045};   // window size in TanY vs iteration - for TMVA
    vector<Double_t> dTanY3B0 = {0.01, 0.02, 0.02, 0.02};       // window size in TanY vs iteration - for zero field
    vector<Double_t> dTanXzero = {0.01, 0.02, 0.02, 0.02};      // window size in TanX vs iteration - for zero field
    // static Double_t dVarX[20] =   {0.004, 0.01, 0.006, 0.02}; // window size in VarX vs iteration - coherent with
    // ptCut AZ-120623 static Double_t ptCut[20] =   {1.0,  0.3,  0.7,  0.1}; // min pxz
    vector<Double_t> ptCut = {1.0, 0.7, 0.3, 0.1};   // min pxz - 120623
    vector<Double_t> ptCutB0 = {1., 1., 1., 1.};     // AZ-010123: min pxz - zero field
    // static Double_t ptCutB0[20] =   {0.5,  0.5, 0.5, 0.5 }; //AZ-010123: min pxz - zero field
    vector<Double_t> curvSta = {
        12.0, 12.0, 12.0, 7.5, 5.5, 4.5,
        3.6,  3.0,  3.0,  3.0, 3.0};   // AZ-310722 - max curvature in stations (for tracks starting in Si)
    vector<Double_t> curvSta41 = {
        12.0, 12.0, 12.0, 12.0, 7.5, 5.5, 4.5,
        3.6,  3.0,  3.0,  3.0,  3.0};   // AZ-190625 - max curvature in stations (for tracks starting in Si) 4Si+VSP
    // static Double_t curvSta4[20] = {12.0, 12.0, 12.0, 10.0, 7.5,  5.5,  4.5,   3.6,   3.0,  3.0,  3.0}; //AZ-051022 -
    // max curvature in stations (for tracks starting in Si) - 4 Si stations
    vector<Double_t> phiXZ = {0.76, 0.69, 0.69, 0.95, 0.74, 0.60, 0.51, 0.44, 0.39, 0.35};   // Max phi_XZ in stations
    // static Double_t phiXZ4[20] =   {0.76, 0.69, 0.69, 0.66, 0.95, 0.74, 0.60,  0.51,  0.44,  0.39,  0.35}; //Max
    // phi_XZ in stations - 4 Si
    vector<Double_t> phiXZ4 = {0.60, 0.63, 0.61, 0.62, 0.95, 0.74,
                               0.60, 0.51, 0.44, 0.39, 0.35};   // AZ-051122 - 4 Si 6_10_14_18
    vector<Double_t> phiXZ41 = {1.00, 0.60, 0.63, 0.61, 0.62, 0.95,
                                0.74, 0.60, 0.51, 0.44, 0.39, 0.35};   // AZ-190625 - 4 Si 6_10_14_18 + VSP
    vector<Double_t> tanXmax = {0.95, 0.83, 0.83, 1.40, 0.91, 0.68, 0.56, 0.47, 0.41, 0.37};   // Max tan_X in stations
    // static Double_t tanXmax4[20] = {0.95, 0.83, 0.83, 0.77, 1.40, 0.91, 0.68,  0.56,  0.47,  0.41,  0.37}; //Max
    // tan_X in stations - 4 Si
    vector<Double_t> tanXmax4 = {0.68, 0.73, 0.69, 0.71, 1.40, 0.91,
                                 0.68, 0.56, 0.47, 0.41, 0.37};   // AZ-051122 - 4 Si 6_10_14_18
    vector<Double_t> tanXmax41 = {1.50, 0.68, 0.73, 0.69, 0.71, 1.40,
                                  0.91, 0.68, 0.56, 0.47, 0.41, 0.37};   // AZ-190625 - 4 Si 6_10_14_18 + VSP
    vector<Double_t> zMean = {18.5, 27.2, 35.8, 60.9, 91.7, 123.4, 154.0, 185.8, 216.5, 248.0};   // Mean Z of stations
    // static Double_t zMean4[20] =   {18.5, 27.2, 35.8, 44.4, 60.9, 91.7, 123.4, 154.0, 185.8, 216.5, 248.0}; //Mean Z
    // of stations - 4 Si AZ-301222 static Double_t zMean4[20] =   {14.9, 23.5, 32.1, 40.8, 60.9, 91.7, 123.4, 154.0,
    // 185.8, 216.5, 248.0}; //AZ-051122 - 4Si 6_10_14_18
    vector<Double_t> zMean4 = {17.6,  26.4,  35.1,  43.7,  63.6, 94.2,
                               126.5, 156.8, 189.0, 219.2, 251.4};   // AZ-301222 - 4Si 6_10_14_18 - beam data
    vector<Double_t> zMean41 = {7.6,  17.6,  26.4,  35.1,  43.7,  63.6,
                                94.2, 126.5, 156.8, 189.0, 219.2, 251.4};   // AZ-190625 - 4Si 6_10_14_18 + VSP
    //*/
 
    FairField* magField = FairRunAna::Instance()->GetField();
    fBy = TMath::Abs(magField->GetBy(0.0, 0.0, 135.0)) / 10;   // max. field in Tesla
 
    fNhitsMin = nHitsMin;
    fdTanX = (fUseTMVA) ? dTanXTMVA : dTanX;
    fdTanY = (fUseTMVA) ? dTanYTMVA : dTanY;
    fdTanY3 = (fUseTMVA) ? dTanY3TMVA : dTanY3;
    fPTcut = ptCut;
    // AZ-311222
    if (fBy < 0.2) {
        fdTanX = dTanXB0;
        fdTanY = dTanYB0;
        fdTanY3 = dTanY3B0;
        fPTcut = ptCutB0;
    }   // AZ-311222
    fdTanXB0 = dTanXzero;   // AZ-010123
    fPhiXZ = (fNSi == 4) ? phiXZ4 : phiXZ;
    if (fNSi == 5)
        fPhiXZ = phiXZ41;   // AZ-190625 - 4Si+VSP
    fTanXmax = (fNSi == 4) ? tanXmax4 : tanXmax;
    if (fNSi == 5)
        fTanXmax = tanXmax41;   // AZ-190625 - 4Si+VSP
    fZmean = (fNSi == 4) ? zMean4 : zMean;
    if (fNSi == 5)
        fZmean = zMean41;   // AZ-190625 - 4Si+VSP
    // fVarx = dVarX;
    fCurvSta = curvSta;
    if (fNSi == 5)
        fCurvSta = curvSta41;
 
    // AZ-240626 for (int j = 0; j < 20; ++j)
    for (int j = 0; j < (int)fPTcut.size(); ++j)
        if (fPTcut[j] > 0.001)
            fPTcut[j] = 1 / fPTcut[j];   // AZ-170722 - curvature
 
    // fExact = 1; // for debug // removed DZ 13.10.2021
    // fExactSel = 13; // for debug
 
    if (fExact) {
        // if (1) {
        //  Open windows
        // Double_t scale = 3;
        Double_t scale = 1;
        // AZ-240626 for (Int_t j = 0; j < 20; ++j) {
        for (Int_t j = 0; j < (int)dTanX.size(); ++j) {
            dTanX[j] *= scale;
            dTanY[j] *= scale;
            // dX[j] *= scale;
        }
    }
 
    // AZ-271222 for (int j = 0; j < 20; ++j) fCurvSta[j] *= (0.8 / fBy); //AZ-290822 - max curvature vs field
    if (fBy > 0.1)
        // AZ-240626 for (int j = 0; j < 20; ++j)
        for (int j = 0; j < (int)fCurvSta.size(); ++j)
            fCurvSta[j] *= (0.8 / fBy);
    // AZ-311222 else fBy = 0.1; // zero mag. field
 
    // Read system info
    fTxBinw = 0.01;
    TString fileGraphs = gSystem->ExpandPathName("$VMCWORKDIR");
    fileGraphs += "/parameters/reco/";
    if (fNSi == 4)
        fileGraphs += "fitGraphs_8kG.txt";   // Run 8
    else
        fileGraphs += "fitGraphs_8kG_vsp.txt";   // AZ-220625 - Run 9
    ReadTxInfo(fileGraphs);
 
    if (fMatBudgetFileName != "")
        fitter.ReadMatBudget(fMatBudgetFileName);
 
    if (fUseTMVA == 2)
        InitTMVA2();   // AZ-221022 - for doublets
    if (fUseTMVA)
        InitTMVA3();   // AZ-300922 - for triplets
    if (fTMVAtracks)
        InitTMVAtracks();   // AZ-040525 - for tracks
 
    // Create tracking detector setups------------------------------------------------------
    TString pathConfig = gSystem->Getenv("VMCWORKDIR");
    fGemStationSet = new BmnGemStripStationSet(pathConfig + "/parameters/gem/XMLConfigs/" + fGemConfigFile);
    fSilStationSet = new BmnSiliconStationSet(pathConfig + "/parameters/silicon/XMLConfigs/" + fSilConfigFile);
 
    cout << "-I- BmnStsVectorFinderV9: Intialisation successfull" << endl;
 
    return kSUCCESS;
}
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::ReadTxInfo(TString fileName)
{
    // Read information about Tx windows
 
    ifstream fin;
    fin.open(fileName);
    if (fin.fail()) {
        cout << " !!! Error opening Tx file " << fileName << endl;
        exit(1);
    }
    TString chline;
    TObjArray* tokens = NULL;
    int ista0 = -1, ista2;
    vector<Float_t> params(10);
 
    while (1) {
        chline.ReadLine(fin);
        if (fin.eof())
            break;
        // cout << chline << endl;
        tokens = chline.Tokenize(" ");
        int ntok = tokens->GetEntriesFast();
        // if (chline.Contains("Station")) {
        if (ntok == 4) {
            // Station info
            int j = 1;
            ista0 = ((TObjString*)tokens->UncheckedAt(j++))->String().Atoi();
            int ista1 = ((TObjString*)tokens->UncheckedAt(j++))->String().Atoi();
            ista2 = 100 * ista0 + ista1;
            if (fStatData.find(ista2) != fStatData.end())
                continue;
            // Float_t txmax = TMath::Abs (((TObjString*)tokens->UncheckedAt(j++))->String().Atof());
            fTxStep = ((TObjString*)tokens->UncheckedAt(j))->String().Atof();
            unordered_map<int, vector<Float_t>> mstat;
            fStatData[ista2] = mstat;
            // fTanXmax[ista0] = TMath::Abs(txmax); // average max Tx (for upper and lower parts)
            fTanXmax[ista0] = -9;
        } else {
            // int ibin = ((TObjString*)tokens->UncheckedAt(0))->String().Atoi();
            Float_t binc = ((TObjString*)tokens->UncheckedAt(0))->String().Atof();
            int ilr = ((TObjString*)tokens->UncheckedAt(1))->String().Atoi();
            Float_t binc1 = binc;
            if (binc1 > 90.0)
                binc1 -= 100.0;
            if (fTanXmax[ista0] < -8) {
                fTanXmax[ista0] = TMath::Abs(binc1 - fTxStep / 2);
                // cout << ista0 << " " << fTanXmax[ista0] << endl;
            }
            int ibin = FindBin(ista0, fTxStep, binc1);
            if (binc > 90.0)
                ibin += 1000;
            for (int j = 0; j < 5; ++j)
                params[ilr * 5 + j] = ((TObjString*)tokens->UncheckedAt(j + 2))->String().Atof();
            // if (ilr > 0) fStatData[ista2].insert(pair<int,vector<Float_t> >(ibin,params));
            if (ilr > 0)
                fStatData[ista2][ibin] = params;
        }
        if (tokens) {
            tokens->Delete();
            delete tokens;
            tokens = NULL;
        }
    }
    fin.close();
    cout << "-I- BmnStsVectorFinderV9: Station info has been read from " << fileName << endl;
}
 
// -----   Public method Exec   --------------------------------------------
void BmnStsVectorFinderV9::Exec(Option_t* opt)
{
 
    TStopwatch sw;
    sw.Start();
 
    //*
    // Reset output array
    if (!fHitArray)
        Fatal("Exec", "No StsHit array");
 
    // AZ-230626
    if (fKillHitsTxP) {
        // Kill hits in momentum space in MC
        KillHitsTxP();
    }
 
    fVectorArray->Delete();
    fMap2Out.clear();   // AZ-210922
    fMap3Out.clear();   // AZ-260922
 
    for (Int_t j = 0; j < 19; ++j) {
        fCandSet[j].clear();
        // fTripleCodes[j].clear();
        fCandCodes[j].clear();
    }
 
    // AZ-100722 Int_t nTracks = fTrackArray->GetEntriesFast();
    // AZ-100722 if (nTracks == 0) return;
    Int_t nTracks = 0;
    if (fTrackArray) {
        nTracks = fTrackArray->GetEntriesFast();
        if (nTracks == 0)
            return;
    }
 
    const Int_t nPass = 4;   // 4; //5; //8;
    Int_t nsta = CbmKF::Instance()->GetNStsStations();
    for (Int_t j = 0; j < nPass; ++j)
        fNskips[j] = nsta - fNhitsMin[j];
 
    Int_t nHits0 = fHitArray->GetEntriesFast(), idmaxP = 0;
    // cout << " !!! nHits0 " << nHits0 << endl;   // AZ-230626
 
    // Fill cluster-to-hit maps
    fHit2id.clear();
    fClusMaps[0].clear();
    fClusMaps[1].clear();
    fmapHits.clear();
 
    for (Int_t ihit = 0; ihit < nHits0; ++ihit) {
        CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ihit);
        if (hit->GetUniqueID() != 999)
            hit->SetUniqueID(0);   // AZ-021125 - reset hit status
        if (hit->GetZ() > 400.0)
            continue;   // AZ-290323 - exclude small GEM downstream
        Int_t iclusF = hit->GetDigi(0);
        Int_t iclusB = hit->GetDigi(1);
        fClusMaps[0].insert(pair<Int_t, Int_t>(iclusF, ihit));
        fClusMaps[1].insert(pair<Int_t, Int_t>(iclusB, ihit));
        Int_t ista = hit->GetStationNr() - 0;
        fNsta = TMath::Max(fNsta, ista);
        // Code below for debugging
        if (fVerbose > 0) {
            set<Int_t> idset = GetHitId(hit, idmaxP);
            for (set<Int_t>::iterator sit = idset.begin(); sit != idset.end(); ++sit)
                fHit2id.insert(pair<Int_t, Int_t>(ihit, *sit));
        }
    }
 
    Int_t nHitsOut = 0;
    discarded = 0;
    // AZ-080623 BuildTrackCand();
    FillHitInfo();   // AZ-080623
 
    for (Int_t ipass = 0; ipass < nPass; ++ipass) {
        // for (Int_t ipass = nPass-1; ipass < nPass; ++ipass) {
        fPass = ipass;
        // Int_t minHits = 5; //4; // minimum number of hits on tracks to accept
        // Int_t minHits = (ipass == 0) ? 15 : 5; // minimum number of hits on tracks to accept
        Int_t minHits = (ipass == 0) ? 15 : -fNhitsMin[ipass - 1];   // minimum number of hits on tracks to accept
        TClonesArray* trArray = (ipass == 0) ? fTrackArray : fVectorArray;
        nHitsOut += ExcludeHits(minHits, trArray);
 
        if (fVerbose > 0) {
            cout << "-I- BmnStsVectorFinderV9: start - " << nHits0 << ", end - " << nHits0 - nHitsOut << endl;
        }
 
        Int_t ntr0 = fVectorArray->GetEntriesFast();
        // BuildTrackCand();
        /*
        //ExtendTracks(fNsta-1);
        FitTracks();
        RemoveFakes();
        RemoveDoubles();
        cout << "\n ***** Pass " << ipass << ": Number of found tracks = "
         << fVectorArray->GetEntriesFast() - ntr0 << " " << fVectorArray->GetEntriesFast() << "\n" << endl;
        //if (fNhitsMin[ipass] <= 4) RemoveFakes();
        */
        BuildDoublets();
        BuildTriplets();
        fTracks.clear();   //
        BuildTracks();
        FitTracks();
        RemoveDoubles();
        if (fVerbose > 0) {
            // if (1) {
            cout << "\n ***** Pass " << ipass << ": Number of found tracks = " << fVectorArray->GetEntriesFast() - ntr0
                 << " " << fVectorArray->GetEntriesFast() << "\n"
                 << endl;
        }
        // GoToTarget(); //AZ-240625
    }
 
    if (fVerbose > 0)
        cout << "discarded " << discarded << " track candidates" << endl;
 
    // Post-processing - try to exclude fake tracks (with too many shared clusters)
    ExcludeFakes();
    if (fTMVAtracks)
        TMVAOutputTra();   // AZ-040525
 
    // AZ-230626
    // cout << " BBBBB hits " << fHitArray->GetEntriesFast() << " " << fVectorArray->GetEntriesFast() << endl;
 
    sw.Stop();
    workTime += sw.RealTime();
}
 
// -------------------------------------------------------------------------
 
Int_t BmnStsVectorFinderV9::ExcludeHits(Int_t minHits, TClonesArray* trArray)
{
    // Exclude hits used for tracking
 
    // AZ-100722 Int_t nTracks = trArray->GetEntriesFast();
    Int_t nTracks = 0;   // AZ-100722
    if (trArray)
        nTracks = trArray->GetEntriesFast();   // AZ-100722
 
    if (nTracks == 0)
        return 0;
    Int_t nHitsOut = 0;
    // AZ-130722 multimap<Int_t,Int_t>::iterator mit;
    // AZ-130722 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
    unordered_multimap<Int_t, Int_t>::iterator mit;
    pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator> ret;   // AZ-130722
 
    for (Int_t itra = 0; itra < nTracks; ++itra) {
        CbmStsTrack* track = (CbmStsTrack*)trArray->UncheckedAt(itra);
        if (track->GetNStsHits() < minHits) {
            trArray->Remove(track);
            continue;
        }
 
        Int_t nHitsTr = track->GetNStsHits();
 
        for (Int_t ihit = 0; ihit < nHitsTr; ++ihit) {
            int indx = track->GetStsHitIndex(ihit);   // AZ-130722
            // AZ-130722 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(track->GetStsHitIndex(ihit));
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);   // AZ-130722
            if (!hit)
                continue;
 
            if (fNhitsMin[fPass] < -5) {
                // if (fNhitsMin[fPass] > 6) { //AZ-180722
                // Exclude all hits created from given clusters
                Int_t iclusF = hit->GetDigi(0);
                ret = fClusMaps[0].equal_range(iclusF);
                for (mit = ret.first; mit != ret.second; ++mit) {
                    CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
                    if (hit1->GetUniqueID() == 0)
                        ++nHitsOut;         // excluded hit - used in track
                    hit1->SetUniqueID(1);   // exclude hit - used in track
                    // AZ-130722 fmapHits[ihit].used = 1;
                    fmapHits[mit->second].used = 1;   // AZ-130722
                }
                Int_t iclusB = hit->GetDigi(1);
                ret = fClusMaps[1].equal_range(iclusB);
                for (mit = ret.first; mit != ret.second; ++mit) {
                    CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
                    if (hit1->GetUniqueID() == 0)
                        ++nHitsOut;         // excluded hit - used in track
                    hit1->SetUniqueID(1);   // exclude hit - used in track
                    // AZ-130722 fmapHits[ihit].used = 1;
                    fmapHits[mit->second].used = 1;   // AZ-130722
                }
            }
            if (hit->GetUniqueID() == 0)
                ++nHitsOut;
            hit->SetUniqueID(1);
            // AZ-130722 fmapHits[ihit].used = 1;
            fmapHits[indx].used = 1;   // AZ-130722
 
            // AZ-130722 Remove used doublets
            //*
            int ista = hit->GetStationNr() - 1;
            ret = fMap2[ista].equal_range(indx);
            fMap2[ista].erase(ret.first, ret.second);
            // AZ-140722 Remove used triplets
            ret = fMap3[ista].equal_range(indx);
            for (mit = ret.first; mit != ret.second; ++mit) {
                candvec cand = fCandVec3[ista][mit->second];
                string code = MakeCode(cand);
                fMapCode3[ista].erase(code);
            }
            fMap3[ista].erase(ret.first, ret.second);
            //*/
        }
    }
 
    trArray->Compress();
    return nHitsOut;
}
 
// -------------------------------------------------------------------------
 
set<Int_t> BmnStsVectorFinderV9::GetHitId(CbmStsHit* hit, Int_t& idmaxP)
{
    // Get IDs contributing to given hit (for debugging purposes)
 
    // AZ-241222 if (fClusArray[1]) return GetHitIdBmn (hit, idmaxP); // BmnRoot data structures
    if (1)
        return GetHitIdBmn(hit, idmaxP);   // BmnRoot data structures - AZ-241222
 
    set<Int_t> ids;
    // return ids; // FIXME
 
    Int_t iclusF = hit->GetDigi(0);
    Int_t iclusB = hit->GetDigi(1);
    CbmStsCluster* clusF = (CbmStsCluster*)fClusArray[0]->UncheckedAt(iclusF);
    CbmStsCluster* clusB = (CbmStsCluster*)fClusArray[0]->UncheckedAt(iclusB);
 
    map<Int_t, Double_t> indF, indB;
    // set<Int_t> ids;
    Int_t nDigis = clusF->GetNDigis();
    Double_t pmax = 0.0;
    TVector3 mom3;
 
    for (Int_t j = 0; j < nDigis; ++j) {
        CbmStsDigiMatch* digiMatch = (CbmStsDigiMatch*)fDigiMatches->UncheckedAt(clusF->GetDigi(j));
        for (Int_t j1 = 0; j1 < 3; ++j1) {
            Int_t ip = digiMatch->GetRefIndex(j1);
            if (ip < 0)
                break;
            CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(ip);
            p->Momentum(mom3);
            indF[p->GetTrackID()] = mom3.Mag();
        }
    }
 
    nDigis = clusB->GetNDigis();
 
    for (Int_t j = 0; j < nDigis; ++j) {
        CbmStsDigiMatch* digiMatch = (CbmStsDigiMatch*)fDigiMatches->UncheckedAt(clusB->GetDigi(j));
        for (Int_t j1 = 0; j1 < 3; ++j1) {
            Int_t ip = digiMatch->GetRefIndex(j1);
            if (ip < 0)
                break;
            CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(ip);
            // Not needed p->Momentum(mom3);
            indB[p->GetTrackID()] = 1.0;   // mom3.Mag();
        }
    }
 
    idmaxP = -1;
 
    for (map<Int_t, Double_t>::iterator mit = indF.begin(); mit != indF.end(); ++mit) {
        Int_t id = (indB.count(mit->first) == 0) ? -mit->first : mit->first;
        ids.insert(id);   // negative ID for ghost crossings
        if (id >= 0 && mit->second > pmax) {
            pmax = mit->second;
            idmaxP = id;
        }
    }
    // This is to be consistent with CBM hit matching procedure
    // (hit is matched with the closest point)
    if (hit->GetRefIndex() >= 0) {
        CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(hit->GetRefIndex());
        idmaxP = p->GetTrackID();
    } else
        idmaxP = -1;
 
    for (map<Int_t, Double_t>::iterator mit = indB.begin(); mit != indB.end(); ++mit) {
        if (ids.count(mit->first) > 0)
            continue;
        ids.insert(-mit->first);   // negative ID for ghost crossings
    }
 
    return ids;
}
 
// -------------------------------------------------------------------------
 
set<Int_t> BmnStsVectorFinderV9::GetHitIdBmn(CbmStsHit* hit, Int_t& idmaxP)
{
    // Get IDs contributing to given hit (for debugging purposes) - for BmnRoot data structures
 
    set<Int_t> ids;
    // int ind = 0; // Si hit
    // if (hit->GetSystemId() == kGEM) ind = 2; // GEM hit
 
    if (hit->GetRefIndex() == -1 || fStsPoints == nullptr || fSilPoints == nullptr)
        ids.insert(-1);   // fake hit
    else {
        // AZ-150722 CbmStsPoint *p = (CbmStsPoint*) fStsPoints->UncheckedAt(hit->GetRefIndex());
        FairMCPoint* p = nullptr;
        if (hit->GetSystemId() == kGEM)
            p = (FairMCPoint*)fStsPoints->UncheckedAt(hit->GetRefIndex());
        else if (hit->GetSystemId() == kVSP)
            p = (FairMCPoint*)fVspPoints->UncheckedAt(hit->GetRefIndex());
        else
            p = (FairMCPoint*)fSilPoints->UncheckedAt(hit->GetRefIndex());
        ids.insert(p->GetTrackID());
    }
    idmaxP = *ids.begin();   // AZ-290922
    return ids;
}
 
// -------------------------------------------------------------------------
 
// AZ-080623 void BmnStsVectorFinder::BuildTrackCand()
void BmnStsVectorFinderV9::FillHitInfo()   // AZ-080623
{
    // Fill hit information
 
    for (Int_t ist = 0; ist < fNsta; ++ist) {
        // AZ-080623 fmapPhx[ist].clear();
        fmapTx[ist].clear();
        fmapTy[ist].clear();
        fmapX[ist].clear();
        fmapY[ist].clear();
        fSeedVec[ist].clear();
        fCandVec[ist].clear();
        // fCandMap2[ist].clear();
        // fCandMap3[ist].clear();
        fCandVec2[ist].clear();
        fMap2[ist].clear();
        fCandVec3[ist].clear();
        fMap3[ist].clear();
        fMapCode3[ist].clear();
        fTxBins[ist].clear();
    }
    // fmapHits.clear();
 
    // AZ-220725 fTxBinw = 0.01;   // bin width for hits
    Int_t nHits = fHitArray->GetEntriesFast(), idmaxP = 0;
    if (fVerbose > 0)
        cout << "nHits " << nHits << endl;
    TVector3 pos;
 
    for (Int_t ih = 0; ih < nHits; ++ih) {
        CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ih);
        if (hit->GetUniqueID())
            continue;   // used hit
        if (hit->GetZ() > 400.0)
            continue;   // AZ-290323 - exclude small GEM downstream
        // if (hit->GetSystemId() == kGEM) continue; //!!!AZ-100123 - exclude GEMs
        // if (hit->GetSystemId() != kGEM) continue; //!!!AZ-100123 - exclude Si
        // else if (hit->GetSystemId() != kGEM && hit->GetStationNr() == 1) continue; //!!!AZ-110123 - exclude station 1
        // else if (hit->GetSystemId() != kGEM && hit->GetStationNr() == 4) continue; //!!!AZ-110123 - exclude station 1
        // if (hit->GetSystemId() == kGEM && hit->GetStationNr() > 7 && TMath::Abs(hit->GetX()) < 10.0 &&
        // TMath::Abs(hit->GetY()) < 15.0) continue; //!!!AZ-130123 - exclude beam region
        /*if (hit->GetSystemId() == kGEM && hit->GetSectorNr() > 4) continue; //300123 - exclude lower GEMs
        else if (hit->GetSystemId() != kGEM) {
           //AZ-300123 - exclude lower Si
           if (hit->GetStationNr() == 1 && hit->GetSectorNr() > 3) continue;
           if (hit->GetStationNr() == 2 && hit->GetSectorNr() > 5) continue;
           if (hit->GetStationNr() == 3 && hit->GetSectorNr() > 7) continue;
           if (hit->GetStationNr() == 4 && hit->GetSectorNr() > 9) continue;
        }*/
        // AZ-260423 - apply cluster charge correlation cut
        // if (hit->GetSystemId() == kGEM && TMath::Abs(hit->GetSignalDiv()+0.08) > 0.5) continue;
        // else if (hit->GetSystemId() != kGEM && TMath::Abs(hit->GetSignalDiv()-0.11) > 0.5) continue;
        if (fExactSel >= 0) {
            // Accept only hits with given ID - for debug
            set<Int_t> ids = GetHitId(hit, idmaxP);
            if (ids.count(fExactSel) == 0)
                continue;
        }
        if (fExact) {
            // Accept only true hits
            set<Int_t> ids = GetHitId(hit, idmaxP);
            if (idmaxP < 0)
                continue;
        }
        Int_t ista = hit->GetStationNr() - 1;
        //!!! For debug
        // if (ista > 2) continue;
        // if (ista >= 4) continue; //AZ-201222 - Si only
        //!!!
        hit->Position(pos);
        Double_t dx = hit->GetX() - fXyzv[0];
        Double_t dy = hit->GetY() - fXyzv[1];
        Double_t dz = hit->GetZ() - fXyzv[2];
        Double_t by = (fBy < 0.2) ? 0.8 : fBy;   // AZ-271222
        // AZ-250822 if (fPass >= 0) fmapHits[ih] = hitinfo(pos,dx/dz,dy/dz);
        // 271222 if (fPass >= 0) fmapHits[ih] = hitinfo (pos, TMath::ATan2(dx,dz)/fPhiXZ[ista]/fZmean[ista]/fBy,
        // dy/dz/fTanXmax[ista]); //AZ-280822 AZ-010123 if (fPass >= 0) fmapHits[ih] = hitinfo (pos,
        // TMath::ATan2(dx,dz)/fPhiXZ[ista]/fZmean[ista]/by, dy/dz/fTanXmax[ista]); //AZ-271222
        Double_t tanx = dx / dz;
        // Double_t tanx = dx / (fZmean[ista] - fXyzv[2]); //AZ-090623 - w.r.t. nominal vertex (mean station Z!!!)
        if (fPass >= 0)
            fmapHits[ih] = hitinfo(pos, TMath::ATan2(dx, dz) / fPhiXZ[ista] / fZmean[ista] / by,
                                   // AZ-070623 dx/dz/fTanXmax[ista], dy/dz/fTanXmax[ista]); //AZ-010123
                                   // tanx/fTanXmax[ista], dy/dz/fTanXmax[ista]); //AZ-070623
                                   // AZ-230623 tanx, dy/dz/fTanXmax[ista]); //AZ-080623
                                   tanx, dy / dz);   // AZ-230623
        fmapX[ista].insert(pair<Double_t, Int_t>(pos[0], ih));
        fmapY[ista].insert(pair<Double_t, Int_t>(pos[1], ih));
        // AZ-080623 fmapPhx[ista].insert(pair<Double_t,Int_t>(fmapHits[ih].phx,ih));
        fmapTx[ista].insert(pair<Double_t, Int_t>(fmapHits[ih].tx, ih));
        fmapTy[ista].insert(pair<Double_t, Int_t>(fmapHits[ih].ty, ih));
        // AZ-070623 - hit binning in Tx
        int ista2 = ista * 100 + ista + 1;
        int bin = FindBin(ista, fTxBinw, fmapHits[ih].tx);
        if (bin < 0)
            continue;   // AZ-270623
        int binStat = FindBin(ista, fTxStep, BinCenterTx(ista, bin, fTxBinw));
        if (dy > 0)
            binStat += 1000;   // AZ-270623
        if (fStatData[ista2].find(binStat) == fStatData[ista2].end())
            continue;
        // int bin = (tanx + fTanXmax[ista]) / fTxBinw;
        if (dy > 0)
            bin += 1000;
        // int bin = (tanxMean + fTanXmax[ista]) / fTxBinw;
        fTxBins[ista].insert(pair<int, int>(bin, ih));
    }   // for (Int_t ih = 0;
    /*
    //const Int_t stastop[9] = { fNsta-1, fNsta-2, fNsta-4}; // stop station vs iteration
    //const Int_t stastop[9] = { fNsta-3, fNsta-2, fNsta-3}; // stop station vs iteration
    Int_t stastop = fNsta - 1 - fNskips[fPass];
    stastop = 0;
 
    if (fVerbose > 0) cout << "fNskips fPass " << fNskips[fPass] << " " << fPass << " stastop " << stastop << " fNsta-1
    " << fNsta-1 << endl;
 
    for (Int_t ista = fNsta-1; ista >= stastop; --ista) {
      if (fVerbose > 0) cout << "xmapsize " << fmapX[ista].size() << endl;
      for (multimap<Double_t,Int_t>::iterator mit = fmapX[ista].begin(); mit != fmapX[ista].end(); ++mit) {
        candvec aaa;
        //AZ-280323 aaa.nskips = fNsta - 1 - ista;
        aaa.nskips = 0; //AZ-280323
        aaa.momxz = 0.0;
        aaa.stahit[ista] = mit->second;
        aaa.code = "-" + to_string(aaa.stahit[ista]) + "-"; // hit index coded
        //if (fExact) {
 
        set<Int_t> ids = GetHitId(mit->second, idmaxP);
        aaa.idmaxP = idmaxP;
      //}
        //cout << "aaa.z " <<((CbmStsHit*) fHitArray->UncheckedAt(aaa.stahit[ista]))->GetZ() << " ";
        fSeedVec[ista].push_back(aaa);
      } //for (multimap<Double_t,Int_t>::iterator mit
 
      if (fVerbose > 0) {
        Int_t ncand = fSeedVec[ista].size();
        cout << " Vector stat: " << ista << " " << ncand;// << endl;
        //AZ-130722 pair<multimap<Int_t, Int_t>::iterator, multimap<Int_t, Int_t>::iterator> ret;
        pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator> ret; //AZ-130722
 
        for (Int_t j = 0; j < ncand; ++j) {
          Int_t ih = fSeedVec[ista][j].stahit[ista];
          ret = fHit2id.equal_range(ih);
          cout << " (" << ih << "*";
          //AZ-130722 for (multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
          for (unordered_multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) { //AZ-130722
            if (mit != ret.first) cout << ":";
            cout << mit->second;
          }
          cout << ")";
        }
        cout << "\n";
      }
    } // for (Int_t ista = fNsta-1;
    */
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::BuildDoublets()
{
    // Build doublets
 
    int rebuild2 = 1;   // AZ-140722
 
    // Run8
    const unordered_map<string, float> scaleMap4 = {
        {"01", 1.9},  {"02", 2.8}, {"12", 1.0}, {"13", 1.5},   {"23", 0.65}, {"24", 1.4}, {"34", 0.95},
        {"35", 1.5},  {"45", 1.0}, {"46", 1.5}, {"56", 0.7},   {"57", 1.0},  {"67", 0.6}, {"68", 0.7},
        {"78", 0.35}, {"79", 0.5}, {"89", 0.3}, {"810", 0.35}, {"910", 0.25}};
    // Run 9
    const unordered_map<string, float> scaleMap5 = {
        {"01", 11.5}, {"02", 15.0}, {"12", 1.9},  {"13", 2.8},  {"23", 1.0},  {"24", 1.5},   {"34", 0.65},
        {"35", 1.4},  {"45", 0.95}, {"46", 1.5},  {"56", 1.0},  {"57", 1.5},  {"67", 0.7},   {"68", 1.0},
        {"78", 0.6},  {"79", 0.7},  {"89", 0.35}, {"810", 0.5}, {"910", 0.3}, {"911", 0.35}, {"1011", 0.25}};
 
    for (Int_t ist = 0; ist < fNsta; ++ist) {
        fCandVec[ist].clear();
        fCandVec2[ist].clear();
        fMap2[ist].clear();
        fCandVec3[ist].clear();
        fMap3[ist].clear();
        fMapCode3[ist].clear();
        // AZ-140722 - Does not seem to speed up
        /*
        if ( (fPass == 0) || (TMath::Abs(fdTanX[fPass]-fdTanX[fPass-1]) > 0.001) || //AZ-140722
         (TMath::Abs(fdTanY[fPass]-fdTanY[fPass-1]) > 0.001) ) {
           fCandVec2[ist].clear();
           fMap2[ist].clear();
           rebuild2 = 1;
        }
        */
    }
    if (rebuild2 == 0)
        return;   // AZ-140722 - do not rebuild doublets
    const unordered_map<string, float>& scaleMap = (fNSi == 4) ? scaleMap4 : scaleMap5;
 
    Int_t idmaxP = 0;
 
    for (Int_t ista = 0; ista < fNsta - 1; ++ista) {
        int nHits = fTxBins[ista].size();
        // cout << " nHits: " << ista << " " << nHits << endl;
        if (nHits == 0)
            continue;
        string ssta0 = to_string(ista);
 
        for (int istanext = ista + 1; istanext < ista + 3; ++istanext) {
            // for (int istanext = ista + 1; istanext < ista + 2; ++istanext) { // no jumps over station
            if (istanext > fNsta - 1)
                continue;
            if (fmapX[istanext].size() == 0)
                continue;
            Double_t dty = fdTanY[fPass];   //, dphx = fdTanX[fPass];
            // dty *= 5; //AZ-040822
            // dtx *= 2; //AZ-040822
            int nskips = istanext - (ista + 1);
            int ibin0 = -99, ista2 = ista * 100 + istanext;
            set<Int_t> setTx;
            Float_t tx = 0;
            string ssta2 = ssta0 + to_string(istanext);
            dty *= (scaleMap.find(ssta2)->second);
 
            for (unordered_multimap<int, int>::iterator mit = fTxBins[ista].begin(); mit != fTxBins[ista].end(); ++mit)
            {
                // AZ-070623 - Loop over binned hits
                int ih = mit->second;
                if (fmapHits[ih].used)
                    continue;   // used hit
                candvec aaa;
                aaa.nskips = 0;
                aaa.momxz = 0.0;
                aaa.stahit[ista] = mit->second;
                aaa.code = "-" + to_string(aaa.stahit[ista]) + "-";   // hit index coded
                set<Int_t> ids = GetHitId(mit->second, idmaxP);
                aaa.idmaxP = idmaxP;
 
                int ibin = mit->first;
                if (ibin != ibin0) {
                    // Compute dTx constraints for this bin
                    ibin0 = ibin;
                    setTx.clear();
                    tx = BinCenterTx(ista, ibin, fTxBinw);
                    Float_t dtx[2] = {0};
                    DTxWindow(ista2, tx, ibin / 1000, fPTcut[fPass], dtx);
                    // dtx[0] *= 2; dtx[1] *= 2; //
                    if (istanext < fNSi) {
                        dtx[0] *= 2;
                        dtx[1] *= 2;
                    }   // AZ-120623
                    // else { dtx[0] *= 1.5; dtx[1] *= 1.5; }
                    // else { dtx[0] *= 1.5; dtx[1] *= 1.5; } // AZ-260625
                    //  Get hits on the downstream station
                    // multimap<Double_t,Int_t>::iterator mitxb = fmapPhx[istanext].lower_bound(phx-dphx);
                    // multimap<Double_t,Int_t>::iterator mitxe = fmapPhx[istanext].upper_bound(phx+dphx);
                    multimap<Double_t, Int_t>::iterator mitxb = fmapTx[istanext].lower_bound(tx + dtx[0]);
                    multimap<Double_t, Int_t>::iterator mitxe = fmapTx[istanext].upper_bound(tx + dtx[1]);
                    for (multimap<Double_t, Int_t>::iterator mit1 = mitxb; mit1 != mitxe; ++mit1)
                        if (fmapHits[mit1->second].used == 0)
                            setTx.insert(mit1->second);
                }
 
                Double_t ty = fmapHits[ih].ty;
                // AZ-230623 ty *= (fTanXmax[ista] / fTanXmax[istanext]); //AZ-280822
                // Double_t tx = fmapHits[ih].tx;
                // tx *= (fTanXmax[ista] / fTanXmax[istanext]); //AZ-010123
                multimap<Double_t, candvec> mapDoublets;   // AZ-251022 TMVA output -> doublets
 
                // Get hits on the downstream station
                multimap<Double_t, Int_t>::iterator mityb = fmapTy[istanext].lower_bound(ty - dty);
                multimap<Double_t, Int_t>::iterator mitye = fmapTy[istanext].upper_bound(ty + dty);
                // Get hits from the acceptance window
                set<Int_t> setTy, intersect;
                for (multimap<Double_t, Int_t>::iterator mit1 = mityb; mit1 != mitye; ++mit1)
                    if (fmapHits[mit1->second].used == 0)
                        setTy.insert(mit1->second);
                set_intersection(setTx.begin(), setTx.end(), setTy.begin(), setTy.end(),
                                 std::inserter(intersect, intersect.begin()));
                // if (fVerbose > 0) cout << " Intersect: " << ista << " " << ibin << " " << setTx.size() << " "
                //             << setTy.size() << " " << intersect.size() << " " << endl;
                Double_t xp = fmapHits[ih].xyz[0];   // AZ-280722
                Double_t zp = fmapHits[ih].xyz[2];   // AZ-280722
                // int n0 = fMap2[ista].size(); //AZ-241022 - for debug
 
                for (set<Int_t>::iterator sit = intersect.begin(); sit != intersect.end(); ++sit) {
                    CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(*sit);
                    if (hit->GetUniqueID())
                        continue;   // used hit
                    if (fmapHits[*sit].used)
                        continue;
                    if (fBy < 0.2 && TMath::Abs(tx - fmapHits[*sit].tx) > fdTanXB0[fPass])
                        continue;   // AZ-010123 - zero mag. field
                    if (fExact) {
                        // Exact track ID match
                        set<Int_t> ids1 = GetHitId(*sit, idmaxP);
                        // if (ids.count(aaa.idmaxP) == 0) continue;
                        if (idmaxP != aaa.idmaxP)
                            continue;
                    }
                    candvec aaa1 = aaa;
                    aaa1.stahit[istanext] = *sit;                   // second hit of the doublet
                    Double_t dx = xp - fmapHits[*sit].xyz[0];       // AZ-280722
                    Double_t dz = zp - fmapHits[*sit].xyz[2];       // AZ-280722
                    aaa1.lengxz = TMath::Sqrt(dx * dx + dz * dz);   // AZ-280722 - doublet length in XZ-plane
                    aaa1.momxz =
                        Curv3(aaa1, aaa1, aaa1, 0);   // AZ-310722 - pxz-estimate with the origin at the primary vertex
                    if (fBy > 0.2
                        && (TMath::Abs(aaa1.momxz) > fPTcut[fPass] || TMath::Abs(aaa1.momxz) > fCurvSta[ista]))
                        continue;   // AZ-010123
                    else if (fBy < 0.2 && TMath::Abs(aaa1.momxz) / fZmean[1] * fZmean[istanext] > fPTcut[fPass])
                        continue;   // AZ-010123 - zero mag. field
                    // AZ-300922 aaa1.ty = fmapHits[*sit].ty; //AZ-290822
                    // AZ-220623 aaa1.ty = ty - fmapHits[*sit].ty; //AZ-300922
                    aaa1.ty = (fmapHits[ih].xyz[1] - fmapHits[*sit].xyz[1]) / aaa1.lengxz;   // AZ-220623
                    aaa1.tx = tx - fmapHits[*sit].tx;   // AZ-010123 - for zero mag. field
                    // AZ-250822 aaa1.varx = varx;
                    // cout << " Doublet: " << aaa1.idmaxP << " " << aaa1.stahit.begin()->first << endl; // debug output
 
                    if (fOut && fOut % 2 == 0) {   // AZ-210922 for writing out
                        string code = to_string(ih) + "-" + to_string(*sit);
                        if (fMap2Out.find(code) == fMap2Out.end())
                            fMap2Out[code] = fPass;
                    }
 
                    // AZ-221022
                    Double_t tmvaOut = 9.9;
                    if (fUseTMVA == 2)
                        tmvaOut = TMVAOutput2(aaa1);
                    // if (fUseTMVA) outf << tmvaOut << " " << fPass << endl;
                    // if (tmvaOut < 0.02) continue;
                    if (tmvaOut < 0.01)
                        continue;
                    aaa1.nskips = nskips;
 
                    mapDoublets.insert(pair<Double_t, candvec>(-tmvaOut, aaa1));   // AZ-251022
                }   // for (set<Int_t>::iterator sit = intersect.begin();
 
                // AZ-251022
                int n2tot = 0, n2max = (fUseTMVA == 2) ? fNbranches * 2 : 999;
 
                for (auto ait = mapDoublets.begin(); ait != mapDoublets.end(); ++ait) {
                    fCandVec2[ista].push_back(ait->second);
                    fMap2[ista].insert(pair<int, int>(ih, fCandVec2[ista].size() - 1));
                    ++n2tot;
                    if (n2tot >= n2max)
                        break;
                }
                // cout << " *** Doublets : " << ista << " " << ih << " " << fMap2[ista].size()-n0 << endl;
            }   // for (unordered_multimap<int,int>::iterator mit = fTxBins[ista].begin();
        }   // for (int istanext = ista + 1;
 
        Int_t ncand = fCandVec2[ista].size();
        if (fVerbose >= 1) {
            cout << " Doublet stat: " << ista << " " << ncand << endl;
            for (int i2 = 0; i2 < ncand; ++i2)
                PrintHits(fCandVec2[ista][i2]);
        }
    }   // for (Int_t ista = 0;
}
 
// -------------------------------------------------------------------------
 
int BmnStsVectorFinderV9::FindBin(int ista, Float_t binw, Float_t tx)
{
    // Find bin number
 
    return (tx + fTanXmax[ista]) / binw;
}
 
// -------------------------------------------------------------------------
 
Float_t BmnStsVectorFinderV9::BinCenterTx(int ista, int ibin, Float_t binw)
{
    // Get Tx bin center
 
    return -fTanXmax[ista] + binw * (ibin % 1000 + 0.5);
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::DTxWindow(int ista2, Float_t tx, int idu, Float_t curvCut, Float_t* dtx)
{
    // Compute dTx window
 
    // void DTxWindow(Double_t curv)
    // Float_t params[2][4] = { {-0.00995846, -81.1826, -181.066, -94.9851},
    //              {0.0032616, 79.4583, 163.907, -415.319} };
    // Float_t tMinMax[2] = {-0.21, 0.1};
    TF1 ff("ff", "pol3", -1, 1);
    TVector3 xyz[3];
 
    int bintx = FindBin(ista2 / 100, fTxStep, tx) + 1000 * idu;
    // if (fStatData[ista2].find(bintx) == fStatData[ista2].end()) {cout << ista2 << " " << tx << " " << bintx << " " <<
    // fStatData[ista2].size() << " " << endl; exit(0);}
    vector<Float_t>& params = fStatData[ista2][bintx];
    Float_t tMinMax[2] = {params[4], params[9]};
 
    for (int ilr = 0; ilr < 2; ++ilr) {
        int ip = ilr * 5;
        ff.SetParameters(params[ip], params[ip + 1], params[ip + 2], params[ip + 3]);
        Float_t cmax = ff.Eval(tMinMax[ilr]);
        // if (cmax <= curv) { cout << tMinMax[ilr] << endl; continue; }
        if (cmax <= curvCut) {
            dtx[ilr] = tMinMax[ilr];
            continue;
        }
        xyz[1][0] = tMinMax[ilr] * curvCut / cmax;
        xyz[1][1] = ff.Eval(xyz[1][0]);
        Float_t dx = tMinMax[ilr] / 10;
        for (int j = -1; j < 2; j += 2) {
            xyz[j + 1][0] = xyz[1][0] + j * dx;
            xyz[j + 1][1] = ff.Eval(xyz[j + 1][0]);
        }
        // xyz[0].Print();
        // xyz[1].Print();
        // xyz[2].Print();
        //  Local parabolic approximation
        TVector3 cba = Parabola(xyz[0], xyz[1], xyz[2]);
        // cba.Print();
        Float_t d = cba[1] * cba[1] - 4 * cba[2] * (cba[0] - curvCut);
        // cout << d << endl;
        Float_t x1 = (-cba[1] + TMath::Sqrt(d)) / 2 / cba[2];
        Float_t x2 = (-cba[1] - TMath::Sqrt(d)) / 2 / cba[2];
        Float_t dthx = x1;
        if (x1 * x2 < 0) {
            if ((ilr == 0 && dthx > 0) || (ilr == 1 && dthx < 0))
                dthx = x2;
        } else
            dthx = TMath::Sign(TMath::Min(TMath::Abs(x1), TMath::Abs(x2)), x1);
        dtx[ilr] = dthx;
        if (ilr == 0)
            dtx[ilr] = TMath::Min(dtx[ilr], -0.01f);
        else
            dtx[ilr] = TMath::Max(dtx[ilr], 0.01f);
        // cout << x1 << " " << x2 << " " << dthx << endl;
    }
}
 
// -------------------------------------------------------------------------
 
Bool_t BmnStsVectorFinderV9::CheckVarx(Double_t dx, Double_t dz, Double_t tx, Double_t distxz, Double_t& varx)
{
    // AZ-310722 - Check Varx
 
    varx = (dx / dz - tx) / distxz;
    // AZ-250822 if (TMath::Abs(varx) > fVarx[fPass]) return kFALSE;
    return kTRUE;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::BuildTriplets()
{
    // Build triplets
 
    // Int_t idmaxP = 0;
 
    // AZ-130722 pair<multimap<Int_t,int>::iterator,multimap<Int_t,int>::iterator> ret;
    // AZ-130722 multimap<Int_t,int>::iterator mit, mit1;
    pair<unordered_multimap<Int_t, int>::iterator, unordered_multimap<Int_t, int>::iterator> ret;   // AZ-130722
    unordered_multimap<Int_t, int>::iterator mit, mit1;                                             // AZ-130722
    // AZ-140725 Double_t pxzCut = (fUseTMVA) ? 0.3 : 0.2; //AZ-041022 - pxz difference cut
    // pxzCut = 0.3; // AZ-220625
    // if (fBy < 0.2) pxzCut = 1.0; //AZ-261222 - zero magnetic field
    const Float_t pxzScale4[19] = {5.0, 3.8, 3.0, 2.0, 1.2, 1.0, 0.9, 1.1, 1.0};        // AZ-260623 - Run 8
    const Float_t pxzScale5[19] = {7.0, 7.0, 3.8, 3.0, 2.0, 1.2, 1.0, 0.9, 1.1, 1.0};   // AZ-220625 - Run 9
    // const Float_t dtymax9[19] = {0.014, 0.019, 0.020, 0.016, 0.009, 0.007, 0.007, 0.006, 0.006, 0.006, 0.006, 0.006};
    // // AZ-260625 - sigma
    const Float_t dtymax9[19] = {0.070, 0.095, 0.100, 0.080, 0.045, 0.035,
                                 0.035, 0.030, 0.030, 0.030, 0.03,  0.03};   // AZ-260625 - 5*sigma
    // const Float_t beamWidth9[19] = {0.19, 0.12, 0.10, 0.10, 0.11, 0.12, 0.14, 0.17, 0.19, 0.22}; // AZ-260624 - sigma
    // const Float_t beamWidth9[19] = {0.95, 0.60, 0.50, 0.50, 0.55, 0.60, 0.70, 0.85, 0.95, 1.10, 1.1, 1.1}; //
    // AZ-260624 - 5*sigma const Float_t beamWidth9[19] = {0.19, 0.12, 0.10, 0.09, 0.10, 0.10, 0.11, 0.13, 0.15, 0.19};
    // // AZ-260624 - sigma, fit with "w"
    const Float_t beamWidth9[19] = {0.95, 0.60, 0.50, 0.45, 0.50, 0.50,
                                    0.55, 0.65, 0.75, 0.95, 1.0,  1.0};   // AZ-260624 - 5*sigma
 
    // const Float_t *pxzScale = (fNSi == 4) ? pxzScale4 : pxzScale5; // AZ-230625
    const Float_t *pxzScale = pxzScale5, *dtymax = dtymax9, *beamWidth = beamWidth9;   // AZ-270625
    if (fNSi == 4) {                                                                   // AZ-270625
        pxzScale = pxzScale4;
        dtymax = &dtymax9[1];
        beamWidth = &beamWidth9[1];
    }
 
    for (Int_t ista = 0; ista < fNsta - 2; ++ista) {
        discarded = 0;
        Int_t nTra = fCandVec2[ista].size();
        // AZ-280323 Int_t istanext = ista + 1;
        // AZ-280323 if (nTra == 0 || fCandVec2[istanext].size() == 0) continue;
        if (nTra == 0)
            continue;
 
        for (mit = fMap2[ista].begin(); mit != fMap2[ista].end(); ++mit) {
            candvec& aaa = fCandVec2[ista][mit->second];   // first doublet
            Int_t istanext = ista + 1 + aaa.nskips;        // AZ-280323
            if (fCandVec2[istanext].size() == 0)
                continue;   // AZ-280323
            // AZ-300922 Double_t tany = fmapHits[aaa.stahit.begin()->second].ty; //AZ-2908822
            // AZ-300922 tany *= (fTanXmax[ista] / fTanXmax[istanext+1]); //AZ-290822
            // Double_t tany = aaa.ty; //AZ-300922
            Double_t tanx = aaa.tx;   // AZ-010123
 
            ret = fMap2[istanext].equal_range(aaa.stahit.rbegin()->second);
            multimap<Double_t, candvec> mapDoublets;   // AZ-110722 second doublets, attached to one first doublet
            int nbr = 0, newtr = 0, newtr3 = 0;        //, mult2 = 0;
            // Double_t c2max = 999.0;
 
            // AZ-090525 - for TMVA
            /*
            if (fUseTMVA && fBy > 0.2) {
            for (mit1 = ret.first; mit1 != ret.second; ++mit1) {
                candvec &aaa2 = fCandVec2[istanext][mit1->second]; // second doublet
                if (TMath::Abs(tany - aaa2.ty) > fdTanY3[fPass]) continue;
                ++mult2; // doublet multiplicity in triplets
            }
            }
            */
 
            for (mit1 = ret.first; mit1 != ret.second; ++mit1) {
                candvec aaa1;
                // AZ-150722 Int_t newtr = (mit1 == ret.first) ? 1 : 0;
 
                // Int_t nhits = 3;
                aaa1 = aaa;
                candvec& aaa2 = fCandVec2[istanext][mit1->second];   // second doublet
                if (fExact) {
                    // Exact track ID match
                    if (aaa.idmaxP != aaa2.idmaxP)
                        continue;
                }
                // AZ-310722 - Check doublet matching criteria
                // if (TMath::Abs(aaa.varx - aaa2.varx) > 0.002) continue; //AZ-050822
                // AZ-070822 if (TMath::Abs(aaa.ty - aaa2.ty) > 0.01) continue; //AZ-040822
                // if (TMath::Abs(aaa.ty - aaa2.ty) > fdTanY3[fPass]) continue; //AZ-250822
                // cout << TMath::Abs(tany - aaa2.ty) << " " << fdTanY3[fPass] << " "
                //      << TMath::Abs(tanx - aaa2.tx) << " " << fdTanXB0[fPass] << endl;
                // AZ-230623 if (TMath::Abs(tany - aaa2.ty) > fdTanY3[fPass]) continue; //AZ-290822
                if (fBy < 0.2 && TMath::Abs(tanx - aaa2.tx) > fdTanXB0[fPass])
                    continue;   // AZ-010123 - for zero mag. field
                                /* AZ-090525
                                        if (fOut / 3 > 0) {   // AZ-260922 for writing out
                                            string code = to_string(aaa.stahit.begin()->second) + "-" +
                                   to_string(aaa.stahit.rbegin()->second)
                                                          + "-" + to_string(aaa2.stahit.rbegin()->second);
                                            if (fMap3Out.find(code) == fMap3Out.end())
                                                fMap3Out[code] = fPass;
                                        }
                                */
 
                /* AZ-260625 - reldiff is not a good criterion for first stations and high momenta for a wide beam
                Float_t reldiff = TMath::Abs(aaa.momxz-aaa2.momxz) / (TMath::Abs(aaa.momxz)+TMath::Abs(aaa2.momxz));
                if (fBy > 0.2 && reldiff > 0.9) continue; //AZ-260623
                if (fBy > 0.2 && reldiff / pxzScale[ista] > pxzCut) continue; //AZ-260623
                //if (fBy < 0.2 && TMath::Abs(aaa.momxz-aaa2.momxz) / TMath::Sqrt(fZmean[ista]) > pxzCut) continue;
                //AZ-311222
                */
                // AZ-280625 if (TMath::Abs (aaa.ty - aaa2.ty) > dtymax[ista]) continue; // AZ-260625
                Double_t dtym = dtymax[ista];                    // AZ-280625
                int dsta = aaa2.stahit.rbegin()->first - ista;   // AZ-280625
                if (dsta > 2)
                    dtym *= (2.0 / dsta);   // AZ-280625
                if (TMath::Abs(aaa.ty - aaa2.ty) > dtym)
                    continue;   // AZ-260625
 
                // AZ-260625 Use parabolic approximation to find X-position of triplet at Ztarg
                map<int, int>::iterator mitt = aaa.stahit.begin();
                TVector3 points[3];
                points[0] = fmapHits[mitt->second].xyz;
                ++mitt;
                points[1] = fmapHits[mitt->second].xyz;
                points[2] = fmapHits[aaa2.stahit.rbegin()->second].xyz;
                for (int j3 = 0; j3 < 3; ++j3)
                    points[j3] -= fXyzv;
                TVector3 cba = Parabola(points[0], points[1], points[2], 2);
                Double_t dist = TMath::Abs(cba[0]) / TMath::Sqrt(points[0].Z());
                // cout << " ddd " << ista << " " << points[0].Z() << " " << points[1].Z() << " " << points[2].Z() << "
                // " << dist << endl;
                Double_t beamW = beamWidth[ista] * 2;   // AZ-290625
                if (dsta > 2)
                    beamW *= (2.0 / dsta);   // AZ-290625
                if (dist > beamW)
                    continue;
 
                aaa1.stahit[aaa2.stahit.rbegin()->first] = aaa2.stahit.rbegin()->second;   // third hit of the triplet
 
                Double_t tmvaOut = 9.9;   // AZ-300922
                                          /* AZ-200525
                                                  if (fUseTMVA)
                                                      tmvaOut = TMVAOutput3(aaa1, aaa2);   // AZ-300922
                                                  // if (fUseTMVA) outf << tmvaOut << " " << fPass << endl; //AZ-300922
                                                  if (tmvaOut < 0.02)
                                                      continue;   // AZ-300922 - use TMVA MLPBNN method
                                          */
 
                map<Int_t, Int_t>& hitMap = aaa1.stahit;
 
                CbmStsHit* hit = nullptr;
                CbmStsTrack track;
                TArrayI& hits = *track.GetStsHits();
                hits.Set(3);
                Int_t indx = 0;   //, iok = 1;
 
                /// hitMap ---- aaa1.stahit
                for (map<Int_t, Int_t>::iterator mit2 = hitMap.begin(); mit2 != hitMap.end(); ++mit2) {
                    hits[indx++] = mit2->second;
                    hit = (CbmStsHit*)fHitArray->UncheckedAt(mit2->second);
                    // cout << " hit: " << mit2->first << " " << mit2->second << endl;
                    // hit->SetDx(0.08/TMath::Sqrt(12.0));
                    hit->SetDx(0.08 / TMath::Sqrt(12.0) * 1.2);   // AZ-210125 - as in vers.1
                    /// hit->SetDx(0.05); //AZ-110123 - early estimate
                    // hit->SetDx(0.015);
                    hit->SetDy(0.1234);   // AZ-210125 - as in vers.1
                    // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
                    if (hit->GetStationNr() <= fNSi)
                        hit->SetDx(0.02 / TMath::Sqrt(12.0));   // AZ-210125 - as in vers.1
                    // if (hit->GetStationNr() <= fNSi) hit->SetDx(0.01); //AZ-040123 - early estimate
                    /// if (hit->GetStationNr() <= fNSi) hit->SetDx(0.02); //AZ-040123 - early estimate
                    // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits
                    // else hit->SetCovXY(-1.968e-3);
                    // hit->SetDy(0.03); // for test
                    if (indx == 3)
                        break;
                }
 
                // fitter.DoFit(&track);
                // float chi2 = track.GetChi2();
                // AZ-310722 Double_t ty = 0.0;
                Double_t ty = aaa.ty;   // AZ-310722
                // AZ-280722 float chi2 = LinearFit (&track, newtr, ty) / hit->GetDy() / hit->GetDy();
                float chi2 = LinearFit(aaa, aaa2, &track, newtr, ty) / hit->GetDy() / hit->GetDy();   // AZ-280722
                ++newtr;                                                                              // AZ-150722
                // float chi2 = LinearFit (&track, newtr, ty) / 0.09 / 0.09;
                // cout << "3hit track Chi2 " << track.GetChi2() << " " << nhits << endl;
                //  if Chi2 > 10.0 exclude track cand
                //  TODO DZ add index so I know how many candidates were discarded
                // if (aaa1.idmaxP == 78 || aaa1.idmaxP == 222 || aaa1.idmaxP == 127)
                // cout << " Triplet: " << aaa1.idmaxP << " " << aaa1.stahit.begin()->first << " " << chi2 << endl; //
                // debug output
 
                if (chi2 > 10.0) {
                    // if (chi2 > 20.0) { //AZ-010123
                    discarded++;
                    continue;
                }
                // Triplet curvature
                aaa1.momxz = Curv3(aaa, aaa2, aaa1, newtr3);
                newtr3 = 1;
                // AZ-150722 if (TMath::Abs(aaa1.momxz) < 0.2) {
                // Double_t ptcut = (fPass < 2) ? 0.3 : 0.2; //AZ-150722
                // if (TMath::Abs(aaa1.momxz) < ptcut) { //AZ-150722
                if (TMath::Abs(aaa1.momxz) > fPTcut[fPass]) {   // AZ-170722 - curvature
                    // if (TMath::Abs(aaa1.momxz) > fPTcut[fPass] || TMath::Abs(aaa1.momxz) > fCurvSta[ista]) {
                    // //AZ-310722 - curvature
                    discarded++;
                    continue;
                }
 
                if (fOut / 3 > 0) {   // AZ-090525 for writing out
                    string code = to_string(aaa.stahit.begin()->second) + "-" + to_string(aaa.stahit.rbegin()->second)
                                  + "-" + to_string(aaa2.stahit.rbegin()->second);
                    if (fMap3Out.find(code) == fMap3Out.end())
                        fMap3Out[code] = fPass;
                }
 
                aaa1.ty = ty;
                /*AZ-110722
                fCandVec3[ista].push_back(aaa1);
                fMap3[ista].insert(pair<Int_t,int>(mit->first,fCandVec3[ista].size()-1));
                string code = "-";
                for (map<Int_t,Int_t>::iterator mitr = hitMap.begin(); mitr != hitMap.end(); ++mitr)
                  code += (to_string(mitr->second) + "-");
                fMapCode3[ista][code] = fCandVec3[ista].size() - 1;
                */
                /*AZ-200525
                        if (fUseTMVA) {
                            // AZ-300922  Keep fNbranches with the highest NN response
                            if (-tmvaOut > c2max && nbr >= fNbranches)
                                continue;
                        } else if (chi2 > c2max && nbr >= fNbranches)
                            continue;
                */
                aaa1.chi2 = chi2;                               // AZ-140722
                Double_t qual = (fUseTMVA) ? -tmvaOut : chi2;   // AZ-300922
                aaa1.nskips = aaa.nskips + aaa2.nskips;         // AZ-280323
                // AZ-300922 mapDoublets.insert(pair<Double_t,candvec>(chi2,aaa1)); //AZ-110722
                mapDoublets.insert(pair<Double_t, candvec>(qual, aaa1));   // AZ-300922
                // mapDoublets.insert(pair<Double_t,candvec>(TMath::Abs(aaa1.momxz),aaa1)); //AZ-280722
                nbr = mapDoublets.size();
                // AZ-310722 if (nbr > fNbranches) mapDoublets.erase(mapDoublets.rbegin()->first);
                /*AZ-200525
                        if (nbr > fNbranches) {
                            multimap<Double_t, candvec>::iterator mit0 = mapDoublets.end();
                            --mit0;
                            mapDoublets.erase(mit0);
                        }
                        c2max = mapDoublets.rbegin()->first;
                */
            }   // for (mit1 = ret.first;
 
            // AZ-200525
            if (fUseTMVA && nbr > fNbranches) {
                multimap<Double_t, candvec> mapTmp;
                for (multimap<Double_t, candvec>::iterator it = mapDoublets.begin(); it != mapDoublets.end(); ++it) {
                    auto aaa1 = it->second;
                    // Double_t tmvaOut = TMVAOutput3 (aaa1, aaa2, mult2); // AZ-090525
                    //  AZ-280525 Double_t tmvaOut = TMVAOutput3 (aaa1, aaa, mult2); // AZ-090525
                    Double_t tmvaOut = TMVAOutput3(aaa1, aaa, nbr);   // AZ-280525
                    if (tmvaOut < -0.2) {
                        ++discarded;
                        continue;
                    }   // AZ-090525 - use TMVA BDTD method
                    mapTmp.insert(pair<Double_t, candvec>(-tmvaOut, aaa1));
                }
                mapDoublets.clear();
                mapDoublets = mapTmp;
            }
            //
 
            // Take only first fNbranches with lowest chi2
            int nok = 0;
            for (multimap<Double_t, candvec>::iterator it = mapDoublets.begin(); it != mapDoublets.end(); ++it) {
                if (nok >= fNbranches)
                    break;
                ++nok;
                fCandVec3[ista].push_back(it->second);
                fMap3[ista].insert(pair<Int_t, int>(mit->first, fCandVec3[ista].size() - 1));
                map<Int_t, Int_t>& hitMap = it->second.stahit;
                string code = "-";
                for (map<Int_t, Int_t>::iterator mitr = hitMap.begin(); mitr != hitMap.end(); ++mitr)
                    code += (to_string(mitr->second) + "-");
                fMapCode3[ista][code] = fCandVec3[ista].size() - 1;
                // cout << ista << " " << code << endl;
            }
        }   // for (mit = fCandMap2[ista].begin();
 
        Int_t ncand = fCandVec3[ista].size();
        if (fVerbose > 0)
            cout << " Triplet stat: " << ista << " " << ncand << " " << discarded << " " << fPass << endl;
    }   // for (Int_t ista = fNsta-1;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::BuildTracks()
{
    // Build tracks starting from triplets
 
    int istaEnd = fNsta - fNhitsMin[fPass] + 1;
    istaEnd = TMath::Min(istaEnd, fNsta - 3);   // do not allow simple triplets on last 3 stations
    // std::multimap<Int_t,int> &candMap = fMap3[ista];
    // vector<candvec> &candVec = fCandVec3[ista];
    // std::multimap<Int_t,int> *candMap = &fMap2[ista];
    // vector<candvec> *candVec = &fCandVec2[ista];
 
    for (Int_t ista = 0; ista < istaEnd; ++ista) {
        // std::multimap<Int_t,int> *candMap = &fMap2[ista];
        // vector<candvec> *candVec = &fCandVec2[ista];
        // AZ-130722 std::multimap<Int_t,int> *candMap = &fMap3[ista];
        std::unordered_multimap<Int_t, int>* candMap = &fMap3[ista];   // AZ-130722
        vector<candvec>* candVec = &fCandVec3[ista];
        Int_t nTra = candMap->size();
        // cout << " BuildTracks: " << ista << " " << nTra << " " << istaEnd << " " << fNsta << endl;
        if (nTra == 0)
            continue;
        // Int_t istaup = ista - 2;
 
        // AZ-130722 for (multimap<Int_t,int>::iterator mit = candMap->begin(); mit != candMap->end(); ++mit) {
        for (unordered_multimap<Int_t, int>::iterator mit = candMap->begin(); mit != candMap->end(); ++mit) {
            candvec& cand = (*candVec)[mit->second];
            if (fPass == 0 && cand.nskips > 0)
                continue;   //!!! AZ-280323 - no jumps over stations at 1'st pass !!!
            map<int, int>::iterator it = cand.stahit.begin();
            int ista0 = it->first;
            ++it;
            int ista1 = it->first;
            if (ista1 - ista0 > 1)
                continue;   // AZ-250423 - isolated first hit is not allowed
            ExtendTrack(cand);
        }
    }
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::ExtendTrack(candvec cand)
{
    // Extend the track by triplets or doublets
 
    const Double_t c2ndfMax = 30.0;   // cut on chi2/NDF
    Double_t pxzCut = 0.3;            // AZ-261222 - pxz difference cut
    // if (fBy < 0.2) pxzCut = 1.0; //AZ-261222 - zero magnetic field
    // const Float_t pxzCuts9pos[19] = {0.043, 0.046, 0.049, 0.060, 0.053, 0.047, 0.050, 0.055, 0.064}; // AZ-290625 -
    // sigma
    const Float_t pxzCuts9pos[19] = {0.215, 0.230, 0.245, 0.300, 0.265,
                                     0.235, 0.250, 0.275, 0.320};   // AZ-270625 - 5*sigma
    // const Float_t pxzCuts9neg[19] = {0.047, 0.043, 0.037, 0.034, 0.032, 0.027, 0.025, 0.024, 0.029}; // AZ-290625 -
    // sigma
    const Float_t pxzCuts9neg[19] = {0.235, 0.215, 0.185, 0.170, 0.160,
                                     0.135, 0.125, 0.120, 0.145};              // AZ-290625 - 5*sigma
    const Float_t* pxzCutsPos = (fNSi == 5) ? pxzCuts9pos : &pxzCuts9pos[1];   // AZ-290625
    const Float_t* pxzCutsNeg = (fNSi == 5) ? pxzCuts9neg : &pxzCuts9neg[1];   // AZ-290625
 
    Int_t nTra = 0, ista = cand.stahit.rbegin()->first;
    // cout << ista << " " << cand.stahit.size() << endl;
    /*
    //!!! Check - step back
    if (ista == fNsta-2) {
      cand.stahit.erase(cand.stahit.rbegin()->first);
      ista = cand.stahit.rbegin()->first;
    }
    */
    // std::multimap<Int_t,int> &candMap = (ista == fNsta-2) ? fMap2[ista] : fMap3[ista]; // doublets or triplets
    // vector<candvec> &candVec = (ista == fNsta-2) ? fCandVec2[ista] : fCandVec3[ista]; // doublets or triplets
    // std::multimap<Int_t,int> *pCandMap = (ista == fNsta-2) ? &fMap2[ista] : &fMap2[ista]; // doublets or triplets
    // vector<candvec> *pCandVec = (ista == fNsta-2) ? &fCandVec2[ista] : &fCandVec2[ista]; // doublets or triplets
    // AZ-130722 std::multimap<Int_t,int> *pCandMap = (ista == fNsta-2) ? &fMap2[ista] : &fMap2[ista]; // doublets or
    // triplets
    std::unordered_multimap<Int_t, int>* pCandMap =
        (ista == fNsta - 2) ? &fMap2[ista] : &fMap2[ista];   // AZ-130722 doublets or triplets
    vector<candvec>* pCandVec = (ista == fNsta - 2) ? &fCandVec2[ista] : &fCandVec2[ista];   // doublets or triplets
    // AZ-130722 multimap<Int_t,int> &candMap = *pCandMap;
    unordered_multimap<Int_t, int>& candMap = *pCandMap;   // AZ-130722
    vector<candvec>& candVec = *pCandVec;
 
    if (ista < fNsta - 1) {
        /*
        if (ista == fNsta-2) {
          // Doublets
          candMap = fMap2[ista];
          candVec = fCandVec2[ista];
        }
        */
        nTra = candMap.count(cand.stahit.rbegin()->second);
    }
 
    // cout << " ExtendTrack " << ista << " " << nTra << endl;
    if (nTra == 0) {
        // No extension found
        if (int(cand.stahit.size()) < fNhitsMin[fPass])
            return;   // too short track
        // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
        // FitTrack(cand);
        string code("-");
        for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
            code += (to_string(it->second) + "-");
        // cout << cand.stahit.rbegin()->first << " " << ista << endl;
        if (fCandCodes[cand.stahit.size()].find(code) == fCandCodes[cand.stahit.size()].end()) {
            fCandCodes[cand.stahit.size()].insert(code);
            // cout << " 1 " << endl;
            // AZ-170122 if (cand.stahit.size() >= 4 && FitTrack(cand) > c2ndfMax) return;
            // AZ-110123 if (cand.stahit.size() >= 4 && cand.stahit.rbegin()->first != ista && FitTrack(cand) >
            // c2ndfMax) return;
            if (int(cand.stahit.size()) >= fNhitsMin[fPass] && cand.stahit.rbegin()->first != ista
                && FitTrack(cand) > c2ndfMax)
                return;   // AZ-110123
            // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
            fCandVec[cand.stahit.begin()->first].push_back(cand);
        }
        return;
    }
 
    // AZ-130722 multimap<Int_t,int>::iterator mit;
    // AZ-130722 pair<multimap<Int_t,int>::iterator,multimap<Int_t,int>::iterator> ret;
    unordered_multimap<Int_t, int>::iterator mit;                                                   // AZ-130722
    pair<unordered_multimap<Int_t, int>::iterator, unordered_multimap<Int_t, int>::iterator> ret;   // AZ-130722
    ret = candMap.equal_range(cand.stahit.rbegin()->second);
    string code = "-";
    int istaMid = 0;
 
    // Use Kalman track constraints around extrapolated position
    /*
    set<int> kalmWind;
    if (nTra > 5)
      if (cand.stahit.size() >= 4) kalmWind = KalmanWindow(cand, candVec[ret.first->second].stahit.begin()->second);
    */
 
    for (mit = ret.first; mit != ret.second; ++mit) {
        // candvec &aaa = candVec[mit->second];
        candvec aaa = candVec[mit->second];
        if (fExact && aaa.idmaxP != cand.idmaxP)
            continue;
        if (fPass == 0 && aaa.nskips > 0)
            continue;   //!!!AZ-289323 - no jumps over station at 1'st pass
 
        // if (nTra > 5 && cand.stahit.size() >= 4 &&
        //  kalmWind.find(aaa.stahit.begin()->second) == kalmWind.end()) continue; // outside Kalman window
 
        // Check if there is a triplet started from the last but one point of the current one (cand)
        map<Int_t, Int_t>::iterator mitr = aaa.stahit.begin();
        ++mitr;
        if (code.size() == 1) {
            map<Int_t, Int_t>::reverse_iterator mit1 = cand.stahit.rbegin();
            ++mit1;
            istaMid = mit1->first;
            code += (to_string(mit1->second) + "-");
            --mit1;
            code += (to_string(mit1->second) + "-");
            code += (to_string(mitr->second) + "-");
        } else {
            // Replace last index
            int pos = code.rfind("-", code.size() - 2);
            code.erase(code.begin() + pos + 1, code.end());
            code += (to_string(mitr->second) + "-");
        }
        // cout << " Check code: " << code << " " << istaMid << endl;
        //  Check extension possibility
        int ok = 1;
        if (fMapCode3[istaMid].find(code) == fMapCode3[istaMid].end())
            ok = 0;   // no triplet found
        //!!! Check triplet slopes and curvatures
        //*
        candvec* bbb = nullptr;
        if (ok) {
            bbb = &fCandVec3[istaMid][fMapCode3[istaMid][code]];
            // Double_t dty = TMath::Abs (cand.ty - bbb.ty);
            // if (dty * 2 / (TMath::Abs(cand.ty) + TMath::Abs(bbb.ty)) > 0.1) continue; // different slopes
            // if (dty > 0.005) continue; // different slopes
            //  Curvatures
            // cout << cand.momxz << " " << bbb.momxz << endl;
            if ((fPass == 0 && bbb->nskips > 0) || (cand.nskips + bbb->nskips > 2))
                ok = 0;   //!!!AZ-280323 - jumps over stations
            if (ok && cand.stahit.size() > 2) {
                Double_t dp = TMath::Abs(cand.momxz - bbb->momxz);
                // if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.3) {
                // AZ-060822 if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.2) {
                // AZ-261222 if (dp / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.3) { //AZ-060822
                //  AZ-270625 if (TMath::Abs(fBy) > 0.2 && dp / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) >
                //  pxzCut) { //AZ-311222
                int istabeg = cand.stahit.begin()->first;                                // AZ-290625
                pxzCut = (cand.momxz > 0) ? pxzCutsPos[istabeg] : pxzCutsNeg[istabeg];   // AZ-290625
                int dsta = bbb->stahit.rbegin()->first - istabeg;                        // AZ-290625
                if (dsta > 3)
                    pxzCut *= (3.0 / dsta);   // AZ-290625
                if (TMath::Abs(fBy) > 0.2
                    && dp / (TMath::Abs(cand.momxz) + TMath::Abs(bbb->momxz)) / TMath::Sqrt(TMath::Abs(1 / cand.momxz))
                           > pxzCut)
                {   // AZ-270625
                    // Different curvatures
                    ok = 0;
                }
            }
        }
        if (ok) {
            if ((fPass == 0 && aaa.nskips > 0) || (cand.nskips + aaa.nskips > 2))
                ok = 0;   //!!!AZ-280323 - jumps over stations
            if (aaa.stahit.size() == 3) {
                // Extension by a triplet
                if (ok) {
                    Double_t dp = TMath::Abs(cand.momxz - aaa.momxz);
                    // if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.3) {
                    // AZ-070822 if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.2) {
                    // AZ-261222 if (dp / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.3) { //AZ-070822
                    //  AZ-270625 if (TMath::Abs(fBy) > 0.2 && dp / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) >
                    //  pxzCut) { //AZ-020123
                    int istabeg = cand.stahit.begin()->first;                                // AZ-290625
                    pxzCut = (cand.momxz > 0) ? pxzCutsPos[istabeg] : pxzCutsNeg[istabeg];   // AZ-290625
                    int dsta = aaa.stahit.rbegin()->first - istabeg;                         // AZ-290625
                    if (dsta > 3)
                        pxzCut *= (3.0 / dsta);   // AZ-290625
                    if (TMath::Abs(fBy) > 0.2
                        && dp / (TMath::Abs(cand.momxz) + TMath::Abs(aaa.momxz))
                                   / TMath::Sqrt(TMath::Abs(1 / cand.momxz))
                               > pxzCut)
                    {   // AZ-270625
                        // different curvatures
                        ok = 0;
                    }
                }
            }
        }
        if (!ok) {
            // No extension found
            if (int(cand.stahit.size()) < fNhitsMin[fPass])
                continue;   // too short track
            // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
            string codeVec("-");
            for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
                codeVec += (to_string(it->second) + "-");
            if (fCandCodes[cand.stahit.size()].find(codeVec) == fCandCodes[cand.stahit.size()].end()) {
                fCandCodes[cand.stahit.size()].insert(codeVec);
                if (cand.stahit.size() >= 4) {
                    // cout << " 2 " << endl;
                    if (cand.stahit.rbegin()->first == ista)
                        continue;   // AZ-171121
                    Double_t chi2 = FitTrack(cand);
                    // cout << " --- Chi2: " << chi2 << endl;
                    if (chi2 > c2ndfMax)
                        continue;
                }
                // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
                fCandVec[cand.stahit.begin()->first].push_back(cand);
            }
            continue;   // AZ-280323
        }
        //*/
        //!!! Extension found
        candvec candext = cand;
        if (cand.stahit.size() < 3)
            candext.momxz = bbb->momxz;   // update momentum
        // candext.ty = bbb.ty;
        //  Add triplet (or doublet)
        for (mitr = aaa.stahit.begin(); mitr != aaa.stahit.end(); ++mitr)
            candext.stahit[mitr->first] = mitr->second;
        candext.nskips += aaa.nskips;   //!!!AZ-280323 - update number of jumps
        // cout << " ids: " << aaa.idmaxP << " " << cand.idmaxP << " " << candext.stahit.size() << " " << candext.nskips
        // << " " << fPass << endl; if (candext.nskips > 2) exit(0); AZ-120722 if (aaa.idmaxP != cand.idmaxP)
        // cand.idmaxP = -1;
        if (aaa.idmaxP != cand.idmaxP)
            candext.idmaxP = -1;   // AZ-120722
        // Fit track to fight with combinatorics
        if (candext.stahit.size() >= 4) {
            Double_t chi2ndf = FitTrack(candext);
            // cout << " -3- " << chi2ndf << endl;
            if (chi2ndf > c2ndfMax)
                continue;
            /* This slows down
            if (chi2ndf > c2ndfMax) {
           string codext = MakeCode(candext);
           fCandCodes[candext.stahit.size()].insert(codext);
           continue;
            }
            */
        }
        ExtendTrack(candext);
    }   // for (mit = ret.first;
}
 
// -------------------------------------------------------------------------
 
Double_t BmnStsVectorFinderV9::FitTrack(candvec& cand)
{
    // Fit track candidate
 
    const Int_t gkChi2Cut = 5.0;   // 10.0;
    // Int_t gkChi2Cut = 10.0; //AZ-010123
 
    CbmStsTrack track;
    string hitcode;
 
    MakeStsTrack(cand, hitcode, track);
    int nhits = track.GetNStsHits();
    if (fCandSet[nhits].find(hitcode) != fCandSet[nhits].end())
        return 0.0;   // track has been saved already
 
    // AZ-160122 fitter.DoFit(&track);
    if (nhits <= 4)
        fitter.DoFit(&track);
    else
        FilterHit(cand, track);
    /*
    else {
      // Debug
      CbmStsTrack track1(track);
      FilterHit (cand, track);
      fitter.DoFit(&track1);
      cout << " debug " << track.GetChi2() << " " << track1.GetChi2() << endl;
    }
    */
    // track.GetParamFirst()->Print();
    // cout << " aaaaaa " << track.GetChi2() << " " << track.GetNStsHits() << endl;
    int ndf = TMath::Max(track.GetNDF(), 1);
    Double_t chi2ndf = track.GetChi2() / ndf;
    cand.param = *(track.GetParamLast());
    cand.chi2 = track.GetChi2();
    if (chi2ndf < gkChi2Cut && nhits >= fNhitsMin[fPass]) {
        // Refit with material budget
        // fitter.DoFit(&track);
        fitter.Fit(&track);   // AZ
        if (track.GetChi2() / ndf > gkChi2Cut) {
            // cout << " yyyyyy " << " " << track.GetChi2() / ndf << endl;
            // fitter.Fit(&track); //AZ - for debugging
            return chi2ndf;
        }
        // Save good tracks
        Double_t qual = nhits;
        qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()), 100.0)) / 101.0;
        fTracks.insert(pair<Double_t, CbmStsTrack>(-qual, track));
        fCandSet[nhits].insert(hitcode);
        if (fVerbose > 0) {
            cout << " FitTrack: " << track.GetChi2() << " ";
            PrintHits(cand);
        }
    }
    // fCandSet[nhits].insert(hitcode); - don't use this !!!
    return chi2ndf;
}
 
// -------------------------------------------------------------------------
//*
Double_t BmnStsVectorFinderV9::FilterHit(candvec& cand, CbmStsTrack& track)
{
    // Filter last hit on track candidate
 
    // const Int_t gkChi2Cut = 10.0;
    // AZ const Int_t gkChi2Cut = 20.0; //AZ-010123
    /*
    CbmStsTrack track;
    string hitcode;
 
    MakeStsTrack(cand, hitcode, track);
    */
    Int_t nhits = cand.stahit.size();
    // CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(track.GetStsHitIndex(nhits-1));
    /*cout << " !!! parameters: " << cand.stahit.size() << " " << cand.param.GetX() << " " << cand.param.GetY()
         << " " << cand.param.GetZ() << " "
         << hit->GetZ() << " " << cand.chi2 << " " << cand.code << endl;*/
    // AZ-050222 fitter.Extrapolate(&track,hit->GetZ(),track.GetParamLast());
    // cout << hit->GetZ() << " " << track.GetParamLast()->GetZ() << endl;
 
    CbmKFTrack kftr;
    kftr.SetStsTrack(track, kFALSE);   // last parameter
    fitter.SetKFHits(kftr, &track);
    Bool_t downstream = kTRUE;
    // Double_t qp0 = 0.0;
    Double_t qp0 = track.GetParamLast()->GetQp();
    // Double_t *cov = kftr.GetCovMatrix();
    // for (Int_t j = 0; j < 15; ++j) cov[j] = 0.0;
    // cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = 10000.0;
    // kftr.GetRefChi2() = 0.0;
    // Double_t chi20 = 0.0, dchi2 = 0.0;
    // Double_t chi20 = kftr.GetRefChi2();
 
    /*CbmKF *KF =*/CbmKF::Instance();
    CbmKFHit* h = NULL;
    // Int_t istold = 0, imax = 0;
    Bool_t err = kFALSE;
 
    // for (Int_t i = 0; i < nhits; i++) {
    h = kftr.GetHit(nhits - 1);
    err = err || h->Filter(kftr, downstream, qp0);
    // Double_t dchi2 = kftr.GetRefChi2() - chi20;
    // cout << err << " " << dchi2 << " " << kftr.GetRefChi2() << " " << cand.chi2+dchi2 << endl;
    kftr.GetTrackParam(*track.GetParamLast());
    track.SetChi2(kftr.GetRefChi2() + err * 1000);
    track.SetNDF(2 * nhits - 5);
    /*
    if (fCandSet[track.GetNStsHits()].find(hitcode) != fCandSet[track.GetNStsHits()].end()) return 0.0; // track has
    been saved already
 
    fitter.DoFit(&track);
    //track.GetParamFirst()->Print();
    //cout << " aaaaaa " << track.GetChi2() << " " << track.GetNStsHits() << endl;
    int ndf = TMath::Max (track.GetNDF(), 1);
    Double_t chi2ndf = track.GetChi2() / ndf;
    if (chi2ndf < gkChi2Cut && track.GetNStsHits() >= fNhitsMin[fPass]) {
      // Save good tracks
      Double_t qual = track.GetNStsHits();
      qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()),100.0)) / 101.0;
      fTracks.insert(pair<Double_t,CbmStsTrack>(-qual,track));
      fCandSet[track.GetNStsHits()].insert(hitcode);
      //cout << " FitTrack: " << hitcode << " " << track.GetChi2() << endl;
    }
    return chi2ndf;
    */
    return 0.0;
}
//*/
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::MakeStsTrack(candvec& cand, string& hitcode, CbmStsTrack& track)
{
    // Fill STS track information
 
    map<Int_t, Int_t>& hitMap = cand.stahit;
    CbmStsHit* hit = nullptr;
    TArrayI& hits = *track.GetStsHits();
    hits.Set(hitMap.size());
    Int_t indx = 0;   //, iok = 1;
    hitcode = "-";
 
    /// hitMap ---- cand.stahit
    for (map<Int_t, Int_t>::iterator mit2 = hitMap.begin(); mit2 != hitMap.end(); ++mit2) {
        hits[indx++] = mit2->second;
        hit = (CbmStsHit*)fHitArray->UncheckedAt(mit2->second);
        // cout << " hit: " << mit2->first << " " << mit2->second << endl;
        hit->SetDx(0.08 / TMath::Sqrt(12.0));   // AZ-210125 - as in vers.1
        /// hit->SetDx(0.05); //AZ-110123 - early estimate
        // hit->SetDx(0.015);
        hit->SetDy(0.1234);   // AZ-210125 - as in vers.1
        // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
        if (hit->GetStationNr() <= fNSi)
            hit->SetDx(0.02 / TMath::Sqrt(12.0));   // AZ-210125 - as in vers.1
        // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits
        // else hit->SetCovXY(-1.968e-3);
        // hit->SetDy(0.03); // for test
        CbmStsStation* stat = CbmStsDigiScheme::Instance()->GetStationByNr(hit->GetStationNr());
        Double_t dx = stat->GetSector(0)->GetDx();
        // if (dx < 0.02) dx *= 2; // scale up for Si
        // AZ-120722 if (dx > 0.01) dx *= 1.2; // for GEMs
        if (dx > 0.02)
            dx *= 1.2;                        // AZ-120722 for GEMs
        hit->SetDx(dx / TMath::Sqrt(12.0));   // AZ-210125 - as in vers.1
        // if (dx < 0.02) hit->SetDx(0.01); //AZ-040123 - early estimate for Si
        /// if (dx < 0.02) hit->SetDx(0.02); //AZ-040123 - early estimate for Si
        hitcode += (to_string(mit2->second) + "-");
    }
    track.SetParamFirst(cand.param);
    track.SetParamLast(cand.param);
    track.SetChi2(cand.chi2);
    cand.code = hitcode;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::ExtendTracks(Int_t ista)
{
    // Extend tracks by going upstream
 
    // const Double_t zsta[19] = {1.132, 32.750, 64.150, 96.750, 129.300, 160.800, 192.100}; // average hit Z-position -
    // Run 6 !!! const Double_t zsta[19] = {16.21, 22.79, 29.22, 42.56, 67.63, 116.10, 137.80, 164.10, 186.20}; //
    // average hit Z-position - Run 7 !!!
    /*
    const Double_t zsta[19] = {18.55, 27.16, 35.78, 61.05, 91.50, 123.50, 153.90, 185.90, 216.40, 248.20}; // average
    hit Z-position - Run 8 !!! const Double_t pars[19][4] = { {1.30358, 20.8953, 0.857369, 0.56339}, {1.64202, 7.18813,
    0.596002, 0.694601}, {3.10765, 3.21564, 0.291252, 1.08062}, {2.25484, 3.81810, 0.430724, 0.840519},
                   {3.38309, 4.33504, 0.349704, 1.07882}, {3.99252, 3.35587, 0.302402, 1.65050},
                   {10.3321, 4.24742, 0.168995, 1.71113} };
    */
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::FinishEvent()
{
    // Finish event
    /*
    FairRootManager* ioman = FairRootManager::Instance();
    */
    if (fStream2out.is_open()) {
        fStream2out << " Event: " << fEvent << " " << fMap2Out.size() << endl;
        for (unordered_map<string, int>::iterator it = fMap2Out.begin(); it != fMap2Out.end(); ++it)
            fStream2out << it->first << " " << it->second << "\n";
    }
    if (fStream3out.is_open()) {
        fStream3out << " Event: " << fEvent << " " << fMap3Out.size() << endl;
        for (unordered_map<string, int>::iterator it = fMap3Out.begin(); it != fMap3Out.end(); ++it)
            fStream3out << it->first << " " << it->second << "\n";
    }
    ++fEvent;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::Finish()
{
    printf("Work time of BmnStsVectorFinderV9: %4.2f sec.\n", workTime);
    if (fStream2out.is_open())
        fStream2out.close();
    if (fStream3out.is_open())
        fStream3out.close();
}
 
//__________________________________________________________________________
 
// AZ-270625 TVector3 BmnStsVectorFinderV9::Parabola(TVector3& pos0, TVector3& pos1, TVector3& pos2)
TVector3 BmnStsVectorFinderV9::Parabola(TVector3& pos0, TVector3& pos1, TVector3& pos2, int ix)   // AZ-270625
{
    // Get parabolic track approximation from 3 points (X vs Z for ix = 2 or Y vs X for ix = 0)
    // y = a*x^2 + bx + c
 
    Double_t x2[3] = {pos0[2] - fXyzv[2], pos1[2] - fXyzv[2], pos2[2] - fXyzv[2]};
    Double_t y2[3] = {pos0[0] - fXyzv[0], pos1[0] - fXyzv[0], pos2[0] - fXyzv[0]};
    Double_t x0[3] = {pos0[0], pos1[0], pos2[0]};
    Double_t y0[3] = {pos0[1], pos1[1], pos2[1]};
    Double_t *x = x0, *y = y0;
    if (ix == 2) {
        x = x2;
        y = y2;
    }
 
    Double_t denom = (x[0] - x[1]) * (x[0] - x[2]) * (x[1] - x[2]);
    Double_t dy10 = y[1] - y[0];
    Double_t dy02 = y[0] - y[2];
    Double_t dy21 = y[2] - y[1];
    // Double_t a = x[2] * (y[1] - y[0]) + x[1] * (y[0] - y[2]) + x[0] * (y[2] - y[1]);
    Double_t a = x[2] * dy10 + x[1] * dy02 + x[0] * dy21;
    a /= denom;
    // Double_t b = x[0]*x[0] * (y[1] - y[2]) + x[2]*x[2] * (y[0] - y[1]) +
    // x[1]*x[1] * (y[2] - y[0]);
    Double_t b = -x[0] * x[0] * dy21 - x[2] * x[2] * dy10 - x[1] * x[1] * dy02;
    b /= denom;
    Double_t c = x[1] * x[1] * (x[2] * y[0] - x[0] * y[2]) + x[1] * (x[0] * x[0] * y[2] - x[2] * x[2] * y[0])
                 + x[0] * x[2] * (x[2] - x[0]) * y[1];
    c /= denom;
    return TVector3(c, b, a);
    // 𝐴=𝑥3(𝑦2−𝑦1)+𝑥2(𝑦1−𝑦3)+𝑥1(𝑦3−𝑦2) / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
    // 𝐵=𝑥21(𝑦2−𝑦3)+𝑥23(𝑦1−𝑦2)+𝑥22(𝑦3−𝑦1) / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
    // 𝐶=𝑥22(𝑥3𝑦1−𝑥1𝑦3)+𝑥2(𝑥21𝑦3−𝑥23𝑦1)+𝑥1𝑥3(𝑥3−𝑥1)𝑦2 / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::FitTracks()
{
    // Fit tracks
 
    const Double_t gkChi2Cut = 5.0;   // 10.0;
    // const Double_t gkChi2Cut = 10.0; //AZ-010123
    // fTracks.clear();
 
    // for (Int_t ista = 0; ista < fNsta-3; ++ista) {
    for (Int_t ista = 0; ista < fNsta; ++ista) {
        Int_t ntra = fCandVec[ista].size();
        if (fVerbose > 0)
            cout << " FitTracks: " << ista << " " << ntra << endl;
 
        for (Int_t itra = 0; itra < ntra; ++itra) {
            candvec& cand = fCandVec[ista][itra];
            // if (cand.nextind >= 0) continue; // track has an extension already processed
 
            map<Int_t, Int_t>& hitMap = cand.stahit;
            Int_t nhits = hitMap.size();
            // if (nhits < 3) continue; // too short track
            if (nhits < fNhitsMin[fPass])
                continue;   // too short track
            CbmStsTrack track;
            TArrayI& hits = *track.GetStsHits();
            hits.Set(nhits);
            Int_t indx = 0, iok = 1;
 
            // Check if the track candidate has been checked already
            string hitcode("-");
            for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit)
                hitcode += (to_string(mit->second) + "-");
            // cout << itra << " " << hitcode << endl;
 
            /* AZ-140725
            Int_t same = 0; //, sta0 = hitMap.begin()->first;
 
            //same = fCandSet[sta0].count(hitcode);
            same = fCandSet[nhits].count(hitcode);
            //if (same) cout << " same: " << sta0 << " " << same << " " << endl;
            */
 
            /////AZ if (same) continue; // used combination - skip it ???
            // fCandSet[sta0].insert(hitcode);
            fCandSet[nhits].insert(hitcode);
 
            for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
                hits[indx++] = mit->second;
                // Debug
                CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
                // cout << " hit: " << mit->first << " " << hit->GetX() << " " << hit->GetY() << " " << hit->GetZ() << "
                // "
                //    << hit->GetDx() << " " << hit->GetDy() << endl;
                // cout << " hit: " << mit->first << " " << mit->second << endl;
                // hit->SetDx(0.08/TMath::Sqrt(12.0));
                /// hit->SetDx(0.05); //AZ-110123 - early estimate
                // hit->SetDx(0.015);
                /// hit->SetDy(0.1234);
                // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
                // if (hit->GetStationNr() <= fNSi) hit->SetDx(0.02/TMath::Sqrt(12.0)); //AZ-030123
                // if (hit->GetStationNr() <= 3) hit->SetDx(0.0025);
                // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits (CbmFitter)
                // else hit->SetCovXY(-1.968e-3);
                // hit->SetDy(0.03); // for test
                CbmStsStation* stat = CbmStsDigiScheme::Instance()->GetStationByNr(hit->GetStationNr());
                Double_t dx = stat->GetSector(0)->GetDx();
                // if (dx < 0.02) dx *= 2; // scale up for Si
                // AZ-120722 if (dx > 0.01) dx *= 1.2; // scale up for GEMs
                if (dx > 0.02)
                    dx *= 1.2;                        // AZ-120722 scale up for GEMs
                hit->SetDx(dx / TMath::Sqrt(12.0));   // AZ-210125 - as in vers.1
                // if (dx < 0.02) hit->SetDx(0.01); //AZ-040123 - early estimate for Si
                /// if (dx < 0.02) hit->SetDx(0.02); //AZ-040123 - early estimate for Si
                /*
                if (fMatBudgetFileName != "" && hit->GetDz() < 1.e-7) {
                  // Store material budget for the hit (along the beam)
                  Double_t zshift = (dx < 0.02) ? 0.5 : 1.5; // Si or GEM
                  TProfile2D *prof = fMatHistos.lower_bound(hit->GetZ()-zshift)->second;
                  int ix = prof->FindBin(hit->GetX());
                  int iy = prof->FindBin(hit->GetY());
                  Double_t radThick = prof->GetBinContent(ix,iy);
                  hit->SetDz(radThick/100);
                }
                */
            }
 
            // cout << " before " << track.GetParamFirst()->GetZ() << endl;
            // fitter.DoFit(&track);
            fitter.Fit(&track);   // AZ
            // cout << " after " << track.GetParamFirst()->GetZ() << " " << track.GetNStsHits() << " " <<
            // track.GetChi2() << endl;
            /*
            CbmKFTrack kftr = CbmKFTrack(track);
            fitter.SetKFHits(kftr, &track);
            kftr.Fit(kTRUE); // downstream
            kftr.Fit(kFALSE); // upstream
            */
            if (fVerbose > 0) {
                cout << " aaaaaaa " << endl;
                // track.GetParamFirst()->Print();
                cout << track.GetChi2() << " " << nhits << " " << endl;
            }
            // AZ-111221 if (track.GetChi2() / track.GetNDF() > gkChi2Cut) {
            if (track.GetChi2() / track.GetNDF() > gkChi2Cut && track.GetNStsHits() > 3) {
                /*
                // Too high chi2/NDF - try to remove one outlier with the largest dChi2
                // Based on the code from CbmKFTrackInterface::Fit()
                CbmKFTrack kftr;
                //kftr.SetStsTrack(track, kFALSE); // last parameter
                //Bool_t downstream = kFALSE;
                kftr.SetStsTrack(track, kTRUE); // first parameter
                fitter.SetKFHits(kftr,&track);
                Bool_t downstream = kTRUE;
                Double_t qp0 = 0.0;
                Double_t *cov = kftr.GetCovMatrix();
                for (Int_t j = 0; j < 15; ++j) cov[j] = 0.0;
                cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = 10000.0;
                kftr.GetRefChi2() = 0.0;
                Double_t chi20 = 0.0, dchi2 = 0.0, dchi2max = 0.0;
 
                CbmKF *KF = CbmKF::Instance();
                CbmKFHit *h = NULL;
                Int_t istold = 0, imax = 0;
                Bool_t err = kFALSE;
                set<Int_t> indok;
 
                //for (Int_t i = nhits - 1; i >= 0; i--) {
                 // h = kftr.GetHit( i );
                 // Int_t ist = h->MaterialIndex;
                 // if (i < nhits - 1) {
                  //  for (Int_t j = istold - 1; j > ist; j--)
                 //     err = err || KF->vMaterial[j]->Pass( kftr, downstream, qp0 );
                 // }
 
                for (Int_t i = 0; i < nhits; i++) {
                  h = kftr.GetHit( i );
                  Int_t ist = h->MaterialIndex;
                  if (i > 0) {
                    for (Int_t j = istold + 1; j < ist; j++)
                      err = err || KF->vMaterial[j]->Pass( kftr, downstream, qp0 );
                  }
                  err = err || h->Filter( kftr, downstream, qp0 );
                  dchi2 = kftr.GetRefChi2() - chi20;
                  cout << dchi2 << " " << kftr.GetRefChi2() << endl;
                  //if (dchi2 <= gkChi2Cut) indok.insert(i);
                  //if (dchi2 <= 2*gkChi2Cut) indok.insert(i);
                  indok.insert(i);
                  if (dchi2 > dchi2max) {
                    dchi2max = dchi2;
                    imax = i;
                  }
                  chi20 = kftr.GetRefChi2();
                  istold = ist;
                }
                */
                // Remove last hit and refit
                set<Int_t> indok;
                for (Int_t i = 0; i < nhits; i++)
                    indok.insert(i);
                int imax = nhits - 1;
 
                if (indok.size() > 3) {
                    // if (indok.size() >= fNhitsMin[fPass]) {
                    // if (indok.size()-1 >= fNhitsMin[fPass]) {
                    indok.erase(imax);
                    // Good short track - refit
                    indx = 0;
                    hitcode = "-";
                    for (set<Int_t>::iterator sit = indok.begin(); sit != indok.end(); ++sit) {
                        hits[indx++] = hits[*sit];
                        // hitcode += hits[*sit];
                        hitcode += (to_string(hits[*sit]) + "-");
                    }
                    // if (fCandSet[sta0].count(hitcode)) continue;
                    // fCandSet[sta0].insert(hitcode);
                    // cout << " failed " << hitcode << endl;
                    if (fCandSet[indok.size()].find(hitcode) != fCandSet[indok.size()].end())
                        continue;
                    fCandSet[indok.size()].insert(hitcode);
 
                    hits.Set(indok.size());
                    // fitter.DoFit(&track);
                    fitter.Fit(&track);   // AZ
                    int ndf = (track.GetFlag() == 0) ? track.GetNDF() : 1;
                    // cout << " bbbbbb " << endl;
                    // track.GetParamFirst()->Print();
                    if (fVerbose > 0)
                        cout << " bbbbbb " << track.GetChi2() << " " << hits.GetSize() << " " << track.GetChi2() / ndf
                             << endl;
                    if (track.GetChi2() / track.GetNDF() > gkChi2Cut)
                        iok = 0;
                } else
                    iok = 0;
 
            }   // if (track.GetChi2() / track.GetNDF() > gkChi2Cut)
 
            // Good track - extra printout
            if (iok) {
                Double_t qual = track.GetNStsHits();
                // qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()),100.0)) / 101.0;
                // if (track.GetNStsHits() <= 4) qual += TMath::Min(1/TMath::Abs(track.GetParamFirst()->GetQp()),10.0)
                // / 11.0;
                if (track.GetNStsHits() < 4) {
                    // Vertex constraint
                    /*
                    TVector3 mom3, pos3;
                    track.GetParamFirst()->Momentum(mom3);
                    track.GetParamFirst()->Position(pos3);
                    //Double_t ang = mom3.Angle(pos3);
                    Double_t phi1 = TMath::ATan2 (track.GetParamFirst()->GetY(), track.GetParamFirst()->GetX());
                    Double_t phi2 = TMath::ATan2 (track.GetParamFirst()->GetTy(), track.GetParamFirst()->GetTx());
                    Double_t dphi = phi1 - Proxim(phi1,phi2);
                    //qual += (2.1 - TMath::Min(TMath::Abs(dphi),2.1)) / 2.2;
                    //Double_t www = TMath::Exp(-TMath::Abs(dphi)/2.1);
                    Double_t www = TMath::Min (TMath::Exp(-TMath::Abs(dphi)/0.2), 0.999); // 0.2 - rms
                    if (TMath::Abs(dphi) > 2.1) www = 0.0;
                    */
                    // Linear extrapolation to Zvert
                    Double_t xextr;   // = pos3.X() - pos3.Z() * track.GetParamFirst()->GetTx() + 0.04;
                    Double_t yextr;   // = pos3.Y() - pos3.Z() * track.GetParamFirst()->GetTy();
                    // Extrapolation to Zvert
                    FairTrackParam param = *track.GetParamFirst();
                    fitter.Extrapolate(&param, fXyzv[2], &param);
                    xextr = param.GetX() - fXyzv[0];
                    yextr = param.GetY() - fXyzv[1];
                    Double_t rad =
                        TMath::Sqrt(xextr * xextr / 0.035 / 0.035 + yextr * yextr / 0.2 / 0.2);   // 0.035, 0.2 - sigmas
                    Double_t www = TMath::Min(TMath::Exp(-rad), 0.999);
                    qual += TMath::Exp(-TMath::Abs(track.GetChi2())) * www;
                } else
                    qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()), 100.0)) / 101.0;
                fTracks.insert(pair<Double_t, CbmStsTrack>(-qual, track));
                // AZ-130722 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
                pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator>
                    ret;   // AZ-130722
                nhits = hits.GetSize();
                if (fVerbose > 0) {
                    cout << " Good track: " << endl;
                    // track.GetParamFirst()->Print();
                    cout << track.GetChi2() << " " << nhits << " " << track.GetParamFirst()->GetZ() << endl;
 
                    for (Int_t j = nhits - 1; j >= 0; --j) {
                        ret = fHit2id.equal_range(hits[j]);
                        cout << " (";
                        // AZ-130722 for (multimap<Int_t,Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
                        for (unordered_multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
                            if (mit == ret.first)
                                cout << mit->first << "*";
                            else
                                cout << ":";
                            cout << mit->second;
                        }
                        cout << ")";
                    }
                    cout << "\n";
                }
            }
 
        }   // for (Int_t itra = 0;
    }
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::RemoveDoubles()
{
    // Remove double tracks (keep the ones with better quality)
 
    multimap<Double_t, CbmStsTrack>::iterator mit = fTracks.begin(), mit1;
    Int_t nOK = fTracks.size(), nOut = fVectorArray->GetEntriesFast();
 
    for (; mit != fTracks.end(); ++mit) {
        if (mit->second.GetFlag())
            continue;   // discarded track
        mit1 = mit;
        ++mit1;
 
        for (; mit1 != fTracks.end(); ++mit1) {
            if (mit1->second.GetFlag())
                continue;   // discarded track
            if (!AreTracksDoubles(mit1->second, mit->second))
                continue;
            mit1->second.SetFlag(1);   // double - discard it
            --nOK;
            /*
            // Debug printout
            TArrayI &hits = *mit1->second.GetStsHits();
            Int_t nhits = mit1->second.GetNStsHits();
            pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
            cout << " Double: " << mit->first << " " << mit1->first << endl;
            for (Int_t j = nhits-1; j >= 0; --j) {
          ret = fHit2id.equal_range(hits[j]);
          cout << " (";
          for (multimap<Int_t,Int_t>::iterator mit2 = ret.first; mit2 != ret.second; ++mit2) {
            if (mit2 == ret.first) cout << mit2->first << "*";
            else cout << ":";
            cout << mit2->second;
          }
          cout << ")";
            }
            cout << "\n";
            */
        }
 
        new ((*fVectorArray)[nOut++]) CbmStsTrack(mit->second);
        // cout << " yyyyy " << mit->second.GetParamFirst()->GetZ() << " " << mit->second.GetNStsHits() << endl;
    }
    if (fVerbose > 0)
        cout << " RemoveDoubles: " << fTracks.size() << " " << nOK << endl;
}
 
// -------------------------------------------------------------------------
 
Bool_t BmnStsVectorFinderV9::AreTracksDoubles(CbmStsTrack& tr1, CbmStsTrack& tr2)
{
    /// Searching common hits in 2 tracks to determine doubles
    /// tr1 contains fewer hits than tr2
 
    Int_t limCommonPoint = (tr1.GetNStsHits() + 1) / 2;   // at least so many common hits should be found
    // limCommonPoint = TMath::Min (limCommonPoint, 2); // 2 common hits at max
    limCommonPoint = TMath::Min(limCommonPoint, 1);   // 1 common hit at max
    TArrayI &hits1 = *tr1.GetStsHits(), &hits2 = *tr2.GetStsHits();
    Int_t nh1 = hits1.GetSize(), nh2 = hits2.GetSize(), nHitsCommon = 0, j = nh2 - 1;
    if (nh1 == 3)
        limCommonPoint = 0;
 
    for (Int_t i = nh1 - 1; i >= 0; i--) {
        // for (Int_t i = 0; i <= nh1 - 1; ++i){
        CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(tr1.GetStsHitIndex(i));
 
        for (; j >= 0; j--) {
            // for ( ; j <= nh2 - 1; ++j){
            CbmStsHit* hit2 = (CbmStsHit*)fHitArray->UncheckedAt(tr2.GetStsHitIndex(j));
 
            // Is the hit common for two tracks?
            if (hit1 == hit2) {
                nHitsCommon++;
                // AZ-281121 if (nHitsCommon == limCommonPoint) return kTRUE; // already enough common hits
                if (nHitsCommon > limCommonPoint)
                    return kTRUE;   // already enough common hits
                break;
            }
 
            if (hit2->GetStationNr() < hit1->GetStationNr())
                break;   // already closer to target
        }
 
        // AZ-281121 if ((nh1 - i - 1) + limCommonPoint - nHitsCommon > nh1) return kFALSE; // there'll be not enough
        // common hits already
        if (i + nHitsCommon <= limCommonPoint)
            return kFALSE;   // there'll be not enough common hits already
    }
 
    // AZ-281121 if (nHitsCommon < limCommonPoint) return kFALSE; // not too many common hits
    if (nHitsCommon <= limCommonPoint)
        return kFALSE;   // not too many common hits
 
    return kTRUE;
}
 
// -------------------------------------------------------------------------
/*
void BmnStsVectorFinderV9::RemoveFakes()
{
  // Remove fake tracks (according to some empirical rules)
 
  Int_t ncand = fVectorArray->GetEntriesFast();
 
  for (Int_t j = 0; j < ncand; ++j) {
    CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(j);
    Int_t nhits = tr->GetNStsHits();
    if (1/TMath::Abs(tr->GetParamFirst()->GetQp()) < 0.2) { fVectorArray->Remove(tr); continue; }
    if (nhits < 4) {
      CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(tr->GetStsHitIndex(0));
      if (hit->GetStationNr() <= 3) continue;
      fVectorArray->Remove(tr); // short track starts in GEM
    }
  }
  fVectorArray->Compress();
 
  cout << " Remove fakes: " << ncand << " " << fVectorArray->GetEntriesFast() << endl;
}
*/
// -------------------------------------------------------------------------
void BmnStsVectorFinderV9::RemoveFakes()
{
    // Remove fake tracks (according to some empirical rules)
 
    multimap<Double_t, CbmStsTrack>::iterator mit = fTracks.begin();
    Int_t nOK = fTracks.size(), nOut = nOK;
 
    for (; mit != fTracks.end(); ++mit) {
        if (mit->second.GetFlag())
            continue;   // discarded track
        CbmStsTrack* tr = &mit->second;
        Int_t nhits = tr->GetNStsHits();
        if (1 / TMath::Abs(tr->GetParamFirst()->GetQp()) < 0.2) {
            tr->SetFlag(1);
            --nOut;
            continue;
        }
        if (nhits < 4) {
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(0));
            if (hit->GetStationNr() <= 3)
                continue;
            tr->SetFlag(1);   // short track starts in GEM
            --nOut;
        }
    }
 
    if (fVerbose > 0)
        cout << " Remove fakes: " << nOK << " " << nOut << endl;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::ExcludeFakes()
{
    // Exclude fake tracks (with too many shared clusters)
 
    Int_t ncand = fVectorArray->GetEntriesFast();
 
    map<Int_t, set<Int_t>> mClusTr;
    map<Int_t, Double_t> mWeight;
 
    for (Int_t j = 0; j < ncand; ++j) {
        CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
        Int_t nhits = tr->GetNStsHits();
 
        for (Int_t jh = 0; jh < nhits; ++jh) {
            Int_t indx = tr->GetStsHitIndex(jh);
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
 
            for (Int_t side = 0; side < 2; ++side) {
                Int_t iclus = hit->GetDigi(side);
                if (mClusTr.count(iclus) == 0) {
                    set<Int_t> aaa;
                    mClusTr[iclus] = aaa;
                }
                mClusTr[iclus].insert(j);
            }
        }
 
        Double_t www = TMath::Exp(-tr->GetChi2() / tr->GetNDF());
        // Extrapolation to Zvert
        if (nhits == 3) {
            FairTrackParam param = *tr->GetParamFirst();
            fitter.Extrapolate(&param, fXyzv[2], &param);
            Double_t xextr = param.GetX() - fXyzv[0];
            Double_t yextr = param.GetY() - fXyzv[1];
            Double_t rad =
                TMath::Sqrt(xextr * xextr / 0.035 / 0.035 + yextr * yextr / 0.2 / 0.2);   // 0.035, 0.2 - sigmas
            www *= TMath::Exp(-rad);
        }
        mWeight[j] = www;
    }
 
    while (1) {
        multimap<Double_t, pair<Int_t, Int_t>> mapQual;
        set<Int_t> set2kill, set2kill1;
 
        for (Int_t j = 0; j < ncand; ++j) {
            CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
            if (tr == NULL)
                continue;   // removed track
            Int_t nhits = tr->GetNStsHits();
            Int_t nover = 0;
 
            for (Int_t jh = 0; jh < nhits; ++jh) {
                Int_t indx = tr->GetStsHitIndex(jh);
                CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
                // if (mClusTr[hit->GetDigi(0)].size() > 1 || mClusTr[hit->GetDigi(1)].size() > 1) ++nover;
 
                for (Int_t side = 0; side < 2; ++side) {
                    Int_t iclus = hit->GetDigi(side);
                    if (mClusTr[iclus].size() > 1)
                        ++nover;
                }
            }
            // Double_t quality = nhits * 1000 + (1 - Double_t(nover) / nhits) * 100
            //+ (1 - TMath::Min(tr->GetChi2(),100.0) / 100);
            Double_t quality = nhits * 1000;
            Double_t www = TMath::Exp(-Double_t(nover) / nhits / 2) * mWeight[j];
            mapQual.insert(pair<Double_t, pair<Int_t, Int_t>>(quality + www, pair<Int_t, Int_t>(j, nover)));
            // if ((nover == nhits) || (nhits == 3 && nover == 2)) set2kill.insert(j);
            // if ((nhits <= 4 && nover == nhits) || (nhits == 3 && nover == 2)) set2kill.insert(j);
            // if (nhits <= 4 && Double_t(nover)/nhits/2 > 0.81) set2kill.insert(j);
            if (nhits <= 5 && Double_t(nover) / nhits / 2 > 0.6)
                set2kill.insert(j);
            else if (!fExact && nhits <= 3 && www < 1.e-8)
                set2kill1.insert(j);   // 1e-8 - empirical value
        }
 
        if (set2kill.size())
            set2kill1.clear();   // first kill in set2kill
        if (set2kill.size() || set2kill1.size()) {
            // Remove track and update cluster-to-track info
 
            for (multimap<Double_t, pair<Int_t, Int_t>>::iterator mit = mapQual.begin(); mit != mapQual.end(); ++mit) {
                Int_t itr = mit->second.first;
                CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
                Int_t nhits = tr->GetNStsHits();   //, nover = mit->second.second;
                if (set2kill.count(itr) == 0 && set2kill1.count(itr) == 0)
                    continue;
 
                for (Int_t jh = 0; jh < nhits; ++jh) {
                    Int_t indx = tr->GetStsHitIndex(jh);
                    CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
                    if (mClusTr[hit->GetDigi(0)].count(itr) == 0 || mClusTr[hit->GetDigi(1)].count(itr) == 0)
                    // AZ-300822 { cout << " !!! No track found: " << itr << endl; exit(0); }
                    {
                        if (fVerbose > 0)
                            cout << " !!! No track found: " << itr << endl; /*exit(0);*/
                    }
                    mClusTr[hit->GetDigi(0)].erase(itr);
                    mClusTr[hit->GetDigi(1)].erase(itr);
                }
                fVectorArray->Remove(tr);
                break;
            }
        } else
            break;
    }   // while (1)
 
    fVectorArray->Compress();
 
    // AZ-280623 - kill short tracks (3 or 4 hits without consecutive triplet)
    int nKilled = 0;
    ncand = fVectorArray->GetEntriesFast();
    /* AZ-190125 - does not kill anything
    for (Int_t j = 0; j < ncand; ++j) {
      CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(j);
      Int_t nhits = tr->GetNStsHits();
      // AZ-190125 if (nhits > 4) continue;
      if (nhits > 44) continue; // AZ-190125
      // AZ-190125 int sta0 = -1, ncons = 0;
      int sta0 = -1, ncons = 1; // AZ-190125
 
      for (Int_t jh = 0; jh < nhits; ++jh) {
        Int_t indx = tr->GetStsHitIndex(jh);
        CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(indx);
        int sta = hit->GetStationNr();
        if (jh == 0) sta0 = hit->GetStationNr();
        if (sta - sta0 < 2) {
       //fVectorArray->RemoveAt(j);
       ++ncons;
        } else ncons = 1;
        if (ncons >= 3) break;
        sta0 = sta;
      }
      if (ncons < 3) { fVectorArray->RemoveAt(j); ++nKilled; }
    }
    fVectorArray->Compress();
    cout << " !!! Killed w/out triplet: " << nKilled << endl;
    */
    // AZ-030723 - check short tracks for missing hits (remove if nmissMax hits in the acceptance are missing)
    ncand = fVectorArray->GetEntriesFast();
    nKilled = 0;
 
    for (Int_t j = 0; j < ncand; ++j) {
        CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
        Int_t nhits = tr->GetNStsHits();
        // if (nhits > 4) continue;
        if (nhits > -4)
            continue;   // AZ-200125 - do not kill
        // if (nhits > 44) continue;
        set<int> setStat;
 
        for (int ih = 0; ih < nhits; ++ih) {
            int indx = tr->GetStsHitIndex(ih);
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
            setStat.insert(hit->GetStationNr());
        }
 
        CbmKFTrack kftr = CbmKFTrack(*tr);
        fitter.SetKFHits(kftr, tr);
        FairTrackParam param;
 
        int firstSta = *setStat.begin(), nmisSi = 0, nmisGem = 0, nmiss = 0;
 
        // Go upstream
 
        for (int ista = firstSta - 2; ista >= 0; --ista) {
            nmiss = 0;
            if (ista >= fNSi) {
                BmnGemStripStation* stat = fGemStationSet->GetStation(ista - fNSi);
                kftr.Smooth(stat->GetZPosition());
                kftr.GetTrackParam(param);
                // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
                // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
                // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
 
                for (int imod = 0; imod < stat->GetNModules(); ++imod) {
                    BmnGemStripModule* mod = stat->GetModule(imod);
                    if (mod->IsPointInsideModule(param.GetX(), param.GetY())) {
                        ++nmiss;
                        ++nmisGem;
                        // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
                        break;
                    }   // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
                }
            } else {
                BmnSiliconStation* stat = fSilStationSet->GetStation(ista);
                kftr.Smooth(stat->GetZPosition());
                kftr.GetTrackParam(param);
                // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
                // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
                // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
 
                for (int imod = 0; imod < stat->GetNModules(); ++imod) {
                    BmnSiliconModule* mod = stat->GetModule(imod);
                    // AZ-180125 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
                    if (mod->IsPointInsideModule(param.GetX(), param.GetY(), kFALSE)) {   // AZ-180125
                        ++nmiss;
                        ++nmisSi;
                        // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
                        break;
                    }   // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
                }
            }
            if (nmiss == 0)
                break;   // no need for checking other stations
        }   // for (int ista = firstSta - 2;
 
        // Check downstream stations
        nmiss = 0;
 
        for (int ista = firstSta; ista <= fNsta; ++ista) {
            if (setStat.find(ista + 1) != setStat.end())
                continue;
            nmiss = 0;
            if (ista >= fNSi) {
                BmnGemStripStation* stat = fGemStationSet->GetStation(ista - fNSi);
                kftr.Smooth(stat->GetZPosition());
                kftr.GetTrackParam(param);
                // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
                // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
                // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
 
                for (int imod = 0; imod < stat->GetNModules(); ++imod) {
                    BmnGemStripModule* mod = stat->GetModule(imod);
                    if (mod->IsPointInsideModule(param.GetX(), param.GetY())) {
                        ++nmiss;
                        ++nmisGem;
                        // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
                        break;
                    }   // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
                }
            } else {
                BmnSiliconStation* stat = fSilStationSet->GetStation(ista);
                kftr.Smooth(stat->GetZPosition());
                kftr.GetTrackParam(param);
                // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
                // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
                // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
 
                for (int imod = 0; imod < stat->GetNModules(); ++imod) {
                    BmnSiliconModule* mod = stat->GetModule(imod);
                    // AZ-180125 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
                    if (mod->IsPointInsideModule(param.GetX(), param.GetY(), kFALSE)) {   // AZ-180125
                        ++nmiss;
                        ++nmisSi;
                        // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
                        break;
                    }   // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
                }
            }
            if (nmiss == 0)
                break;   // no need for checking other stations
        }   // for (int ista = firstSta; ista <= fNsta;
 
        if (nmisSi >= 2) {
            fVectorArray->RemoveAt(j);
            ++nKilled;
            continue;
        }
        tr->SetFlag(nmisSi * 10 + nmisGem);
 
    }   // for (Int_t j = 0; j < ncand;
    // cout << " !!! Killed due to missing: " << nKilled << endl;
    fVectorArray->Compress();
 
    if (fVerbose > 0)
        cout << " Exclude fakes: " << ncand << " " << fVectorArray->GetEntriesFast() << endl;
}
 
// -------------------------------------------------------------------------
 
Double_t BmnStsVectorFinderV9::DxVsMom(Int_t ista, candvec& vec)
{
    // Compute momentum correction for parabolic extrapolation
 
    Double_t pars[19][4] = {{0, 0, 0, 0},
                            {0, 0, 0, 0},
                            {8.28698e-06, 13.3515, 0, 0},
                            {0.031189, 4.56933, 0.031189, 4.56933},
                            {0.395468, 3.77954, 0.190848, 3.41027},
                            {0.128467, 4.26652, 0.128467, 4.26652},
                            {0.160979, 4.95095, 0.160979, 4.95095}};
 
    Double_t pxz = vec.momxz;
    Int_t endsta = vec.stahit.rbegin()->first;
    Double_t p0 = pars[ista][0];
    Double_t p1 = pars[ista][1];
    if (endsta > ista + 2) {
        p0 = pars[ista][2];
        p1 = pars[ista][3];
    }
 
    Double_t dx = p0 / TMath::Power(TMath::Log10(10 * pxz), p1);
    return dx;
}
 
// -------------------------------------------------------------------------
 
Double_t BmnStsVectorFinderV9::Proxim(Double_t phi0, Double_t phi)
{
    /// Adjust angle phi to be "around" phi0 - to avoid discontinuity around +- Pi
 
    Double_t dPhi = phi0 - phi;
    if (TMath::Abs(dPhi) > TMath::Pi())
        phi += TMath::Pi() * 2 * TMath::Sign(1., dPhi);
    return phi;
}
 
// -------------------------------------------------------------------------
// AZ-280722 Double_t BmnStsVectorFinderV9::LinearFit(CbmStsTrack *tr, Int_t newtr, Double_t& ty)
Double_t BmnStsVectorFinderV9::LinearFit(candvec& cand,
                                         candvec& cand2,
                                         CbmStsTrack* tr,
                                         Int_t newtr,
                                         Double_t& ty)   // AZ-280722
{
    /// Linear LS fit of 3 points (in R-Y plane): y = a + b * l, l is a track projection length on X-Z plane
 
    // static Double_t x[3], y[3], z[3], xs[3], ys[3], xys[3], x2s[3];
    static Double_t /*x[3],*/ y[3], /*z[3],*/ l[3] = {0}, lens[3] = {0}, ys[3], lys[3], l2s[3];
 
    // AZ-150722 Int_t i0 = (newtr == 1) ? 0 : 2;
    Int_t i0 = (newtr == 0) ? 0 : 2;
 
    for (Int_t i = i0; i < 3; ++i) {
        CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(i));
        // AZ-280722 x[i] = hit->GetX();
        // AZ-280722 z[i] = hit->GetZ();
        if (i) {
            /*AZ-280722
           Double_t dx = x[i] - x[i-1];
           Double_t dz = z[i] - z[i-1];
           l[i] = lens[i] = TMath::Sqrt (dx * dx + dz * dz) + l[i-1];
           */
            if (i == 1)
                l[i] = lens[i] = cand.lengxz;   // AZ-280722
            else
                l[i] = lens[i] = cand2.lengxz + l[i - 1];   // AZ-280722
        }
        y[i] = ys[i] = hit->GetY();
        // cout << " ily " << i << " " << l[i] << " " << y[i] << endl;
    }
 
    for (Int_t i = i0; i < 3; ++i) {
        lys[i] = l[i] * y[i];
        l2s[i] = l[i] * l[i];
        if (i) {
            lens[i] += lens[i - 1];
            ys[i] += ys[i - 1];
            lys[i] += lys[i - 1];
            l2s[i] += l2s[i - 1];
        }
    }
    Double_t b = 3 * lys[2] - lens[2] * ys[2];
    b /= (3 * l2s[2] - lens[2] * lens[2]);
    Double_t a = (ys[2] - b * lens[2]) / 3;
 
    Double_t chi2 = 0.0;
 
    for (int i = 0; i < 3; ++i) {
        Double_t dy = a + b * l[i] - y[i];
        chi2 += (dy * dy);
    }
    // Check
    // TGraph gr(3,l,y);
    // gr.Fit("pol1");
    // cout << " LS fit: " << a << " " << b << " " << chi2 << endl;
    ty = b;
    return chi2;
}
 
// -------------------------------------------------------------------------
 
Double_t BmnStsVectorFinderV9::Curv3(candvec& cand1, candvec& cand2, candvec& cand3, int newtr)
{
    // Compute triplet curvature in ZX plane
 
    static TVector3 points3[3], midPoint;
    int indx = 1;
    static Double_t by = 0.0;
    // static int first = 99;
 
    // if (first != fPass) {
    //  first = fPass;
    //  AZ-270625 if (first) {
    // if (cand3.stahit.size() <= 2 && first) { // AZ-270625 - for doublets
    if (cand3.stahit.size() <= 2) {   // AZ-270625 - for doublets
        // first = 0;
        points3[0].SetXYZ(fXyzv[0], fXyzv[1], fXyzv[2]);
    }
 
    map<int, int>::iterator mit = cand3.stahit.begin();
    if (0) {   // fPass > 1) { //0 - assuming tracks coming from (0,0,0)
        // if (newtr == 0 && cand3.stahit.size() > 2) {
        points3[0] = fmapHits[mit->second].xyz;   // AZ-280722
        points3[0].SetY(0.0);                     // AZ-280722
    }
    // AZ-310722 ++mit;
    //  AZ-270625 if (cand3.stahit.size() > 2) ++mit; //AZ-310722
    if (cand3.stahit.size() > 2) {
        indx = 0;
        // first = 99;
    }   // AZ-270625
 
    for (; mit != cand3.stahit.end(); ++mit) {
        // if (newtr > 0 && indx < 2) { ++indx; continue; }
        points3[indx] = fmapHits[mit->second].xyz;
        if (indx == 1)
            midPoint = points3[indx];
        // AZ-270625 points3[indx] -= points3[0];
        if (indx)
            points3[indx] -= points3[0];   // AZ-270625
        points3[indx++].SetY(0.0);
    }
 
    TVector3 vec21 = points3[2] - points3[1];
    // AZ-280722 Double_t cosAlpha = points3[2] * vec21 / points3[2].Mag() / vec21.Mag();
    // Double_t cosAlpha = points3[1] * vec21 / points3[1].Mag() / vec21.Mag(); //AZ-280722
    Double_t cosAlpha = points3[1] * vec21 / points3[1].Mag() / cand2.lengxz;   // AZ-280722
    // AZ-280722 Double_t rad = points3[1].Mag() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
    Double_t rad = points3[2].Mag() / 2. / TMath::Sin(TMath::ACos(cosAlpha));   // AZ-280722
    // Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
    // Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
    // AZ-310722 Double_t sign = TMath::Sign (1.0, points3[1].Cross(points3[2]).Y());
    Double_t sign = TMath::Sign(3e-4, points3[1].Cross(points3[2]).Y());   // AZ-310722
    if (newtr == 0) {
        FairField* magField = FairRunAna::Instance()->GetField();
        by = magField->GetBy(midPoint[0], midPoint[1], midPoint[2]);
        if (TMath::Abs(by) < 0.1)
            by = -10.0;   // AZ-261222 - zero magnetic field
    }
    // AZ-270625 Double_t pxz = by * rad * sign;
    Double_t pxz = -by * rad * sign;
    // cout << " pt " << pxz << endl;
    // AZ-170722 return pxz;
    return 1 / pxz;   // AZ-170722
}
// -------------------------------------------------------------------------
 
set<int> BmnStsVectorFinderV9::KalmanWindow(candvec& cand, int hitIndx)
{
    // Find hits around extrapolated track
 
    CbmStsTrack track;
    string hitcode;
    MakeStsTrack(cand, hitcode, track);
    fitter.DoFit(&track);
    FairTrackParam param;
    CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(hitIndx);
    Double_t z = hit->GetZ();
    fitter.Extrapolate(&track, z, &param);
    Double_t sigx = 5 * TMath::Sqrt(param.GetCovariance(0, 0));
    Double_t sigy = 5 * TMath::Sqrt(param.GetCovariance(1, 1));
    if (fVerbose > 0)
        cout << " Kalman: " << z << " " << hit->GetStationNr() << " " << sigx << " " << sigy << endl;
 
    // Get hits on the downstream station
    int istanext = hit->GetStationNr() - 1;
    multimap<Double_t, Int_t>::iterator mityb = fmapY[istanext].lower_bound(param.GetY() - sigy);
    multimap<Double_t, Int_t>::iterator mitye = fmapY[istanext].upper_bound(param.GetY() + sigy);
    multimap<Double_t, Int_t>::iterator mitxb = fmapX[istanext].lower_bound(param.GetX() - sigx);
    multimap<Double_t, Int_t>::iterator mitxe = fmapX[istanext].upper_bound(param.GetX() + sigx);
    // Get hits from the acceptance window
    set<Int_t> setX, setY, intersect;
    for (multimap<Double_t, Int_t>::iterator mit = mitxb; mit != mitxe; ++mit)
        if (fmapHits[mit->second].used == 0)
            setX.insert(mit->second);
    for (multimap<Double_t, Int_t>::iterator mit = mityb; mit != mitye; ++mit)
        if (fmapHits[mit->second].used == 0)
            setY.insert(mit->second);
    set_intersection(setX.begin(), setX.end(), setY.begin(), setY.end(), std::inserter(intersect, intersect.begin()));
    return intersect;
}
 
// -------------------------------------------------------------------------
 
std::string BmnStsVectorFinderV9::MakeCode(candvec& cand)
{
    // Make track code
 
    string code("-");
    for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
        code += (to_string(it->second) + "-");
    return code;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::PrintHits(candvec& cand)
{
    // Print hits
 
    for (map<int, int>::iterator mit = cand.stahit.begin(); mit != cand.stahit.end(); ++mit) {
        int ista = mit->first, ih = mit->second;
        CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ih);
        // cout << mit->first << "(" << hit->GetZ() << "):" << mit->second << " - ";
        printf("%d(%5.2f,%5.2f,%5.2f):%d", ista, hit->GetZ(), fmapHits[ih].tx, fmapHits[ih].ty, ih);
        auto ret = fHit2id.equal_range(ih);
        for (unordered_multimap<Int_t, Int_t>::iterator mit1 = ret.first; mit1 != ret.second; ++mit1) {
            if (mit1 == ret.first)
                cout << "*";
            else
                cout << ":";
            cout << mit1->second;
        }
        cout << " - ";
    }
    cout << endl;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::InitTMVA3()
{
    // Initialize TMVA package for triplets - AZ-300922
 
    // This loads the library
    // TMVA::Tools::Instance();
 
    // Default MVA methods to be trained + tested
    std::map<std::string, int> Use;
 
    // Use["MLPBNN"] = 1;   // Recommended ANN with BFGS training method and bayesian regulator
    Use["BDTD"] = 1;   // decorrelation + Adaptive Boost
    Float_t* var = fVarTMVA;
 
    fReadersTMVA3.insert(pair<int, TMVA::Reader*>(-13, nullptr));   // pass3 for negative
    fReadersTMVA3.insert(pair<int, TMVA::Reader*>(-11, nullptr));
    fReadersTMVA3.insert(pair<int, TMVA::Reader*>(-10, nullptr));
    fReadersTMVA3.insert(pair<int, TMVA::Reader*>(10, nullptr));
    fReadersTMVA3.insert(pair<int, TMVA::Reader*>(11, nullptr));
    fReadersTMVA3.insert(pair<int, TMVA::Reader*>(13, nullptr));
 
    // --- Create the Reader objects
 
    for (auto mit = fReadersTMVA3.begin(); mit != fReadersTMVA3.end(); ++mit) {
        TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
        fReadersTMVA3[mit->first] = reader;
 
        // Create a set of variables and declare them to the reader
        // The variable names MUST corresponds in name and type to those given in the weight file(s) used
 
        // FIXME
        int iii = 0;
 
        reader->AddVariable("tanx", &var[iii++]);
        reader->AddVariable("tany", &var[iii++]);
        reader->AddVariable("zb", &var[iii++]);
        reader->AddVariable("dty", &var[iii++]);
        // reader->AddVariable( "1/pxz1",                 &var[iii++] );
        // reader->AddVariable( "log(abs(pxz1))",         &var[iii++] );
        reader->AddVariable("log(abs(pxz2))", &var[iii++]);
        reader->AddVariable("dpxz", &var[iii++]);
        reader->AddVariable("corrq[0]", &var[iii++]);
        reader->AddVariable("corrq[1]", &var[iii++]);
        reader->AddVariable("corrq[2]", &var[iii++]);
        // reader->AddVariable( "log(mult+rndm)",         &var[iii++] );
        reader->AddVariable("log(mult)", &var[iii++]);
 
        // --- Book the MVA methods
 
        TString dir = "datasetTMVA3/weights";
        dir += mit->first;
        dir += "/";
        TString prefix = "TMVAClassification";
 
        // Book method(s)
        for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
            if (it->second) {
                TString methodName = TString(it->first) + TString(" method");
                TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
                reader->BookMVA(methodName, weightfile);
            }
        }
    }   // for (auto mit = fReadersTMVA3.begin();
}
 
// -------------------------------------------------------------------------
 
Float_t BmnStsVectorFinderV9::TMVAOutput3(candvec& aaa1, candvec& aaa2, int mult)
{
    // Compute TMVA output for triplets: aaa1 - triplet, aaa2 - first doublet
 
    Float_t* var = fVarTMVA;
    map<Int_t, Int_t>& hitMap = aaa1.stahit;
    CbmStsHit* hit[3] = {nullptr};
    int iii = 0, indx = 0, inds[3];
    inds[0] = inds[1] = inds[2] = 0;
 
    /// hitMap ---- aaa.stahit
    for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
        inds[indx] = mit->second;
        hit[indx++] = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
    }
 
    Double_t dx = hit[1]->GetX() - hit[0]->GetX();
    Double_t dz = hit[1]->GetZ() - hit[0]->GetZ();
    Double_t lenxz1 = TMath::Sqrt(dx * dx + dz * dz);
    Double_t dx1 = hit[2]->GetX() - hit[1]->GetX();
    Double_t dz1 = hit[2]->GetZ() - hit[1]->GetZ();
    Double_t lenxz2 = TMath::Sqrt(dx1 * dx1 + dz1 * dz1);
    Double_t dty = (hit[2]->GetY() - hit[1]->GetY()) / lenxz2 - (hit[1]->GetY() - hit[0]->GetY()) / lenxz1;
    Double_t pxz1 = aaa1.momxz;   // triplet curvature
    Double_t dpxz = (aaa2.momxz - pxz1) / (TMath::Abs(aaa2.momxz) + TMath::Abs(pxz1));
 
    var[iii++] = hit[0]->GetX() / hit[0]->GetZ();
    var[iii++] = hit[0]->GetY() / hit[0]->GetZ();
    var[iii++] = hit[0]->GetZ();
    var[iii++] = dty;
    // var[iii++] = pxz1;
    var[iii++] = TMath::Log(TMath::Abs(1 / pxz1));
    var[iii++] = dpxz;
    var[iii++] = hit[0]->GetSignalDiv();
    var[iii++] = hit[1]->GetSignalDiv();
    var[iii++] = hit[2]->GetSignalDiv();
    var[iii++] = TMath::Log(mult);
 
    int iread = (fPass == 2) ? 1 : fPass;   // fPass==2 has the same triplet cuts as fPass==1
    iread += 10;
    if (pxz1 < 0)
        iread = -iread;
    Float_t tmvaout = fReadersTMVA3[iread]->EvaluateMVA("BDTD method");
    /*
    if (fOutf) {
    for (int j = 0; j < iii; ++j) *fOutf << var[j] << " ";
    *fOutf << mult << " " << tmvaout << " " << iread << " " << inds[0] << " " << inds[1] << " "
           << inds[2] << " " << fEvNo << endl;
    }
    */
    return tmvaout;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::InitTMVA2()
{
    // Initialize TMVA package for doublets - AZ-221022
 
    // This loads the library
    // TMVA::Tools::Instance();
 
    // Default MVA methods to be trained + tested
    std::map<std::string, int> Use;
 
    Use["MLPBNN"] = 1;   // Recommended ANN with BFGS training method and bayesian regulator
    Use["BDTD"] = 1;     // decorrelation + Adaptive Boost
 
    // --- Create the Reader object
 
    TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
    fReaderTMVA2 = reader;
 
    // Create a set of variables and declare them to the reader
    // The variable names MUST corresponds in name and type to those given in the weight file(s) used
 
    // FIXME
    Float_t* var = fVarTMVA;
    int iii = 0;
    reader->AddVariable("xb", &var[iii++]);
    reader->AddVariable("yb", &var[iii++]);
    reader->AddVariable("zb", &var[iii++]);
    reader->AddVariable("tx", &var[iii++]);
    reader->AddVariable("ty", &var[iii++]);
    // reader->AddVariable( "dphixsc",              &var[iii++] );
    reader->AddVariable("dtanysc", &var[iii++]);
    reader->AddVariable("lenxz", &var[iii++]);
    reader->AddVariable("1/pxz", &var[iii++]);
    reader->AddVariable("corrq[1]-corrq[0]", &var[iii++]);
    reader->AddVariable("corrn[1]-corrn[0]", &var[iii++]);
    reader->AddVariable("pass", &var[iii++]);
 
    // --- Book the MVA methods
 
    TString dir = "dataset4diff/weights2/";
    TString prefix = "TMVAClassification";
 
    // Book method(s)
    for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
        if (it->second) {
            TString methodName = TString(it->first) + TString(" method");
            TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
            reader->BookMVA(methodName, weightfile);
        }
    }
}
 
// -------------------------------------------------------------------------
 
Float_t BmnStsVectorFinderV9::TMVAOutput2(candvec& aaa1)
{
    // Compute TMVA output for doublets
 
    Float_t* var = fVarTMVA;
    map<Int_t, Int_t>& hitMap = aaa1.stahit;
    CbmStsHit *hit = nullptr, *hit1 = nullptr;
    int iii = 0, indx = 0;
 
    /// hitMap ---- aaa.stahit
    for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
        hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
        if (indx == 0) {
            var[iii++] = hit->GetX();
            var[iii++] = hit->GetY();
            var[iii++] = hit->GetZ();
            hit1 = hit;
        } else if (indx == 1) {
            Double_t dx = hit->GetX() - var[0];
            Double_t dy = hit->GetY() - var[1];
            Double_t dz = hit->GetZ() - var[2];
            var[iii++] = dx / dz;
            var[iii++] = dy / dz;
            var[iii++] = aaa1.ty;
            var[iii++] = aaa1.lengxz;
            var[iii++] = aaa1.momxz;
            var[iii++] = hit->GetSignalDiv() - hit1->GetSignalDiv();
            var[iii++] = hit->GetTimeStamp() - hit1->GetTimeStamp();
            var[iii++] = (fPass == 2) ? 1 : fPass;   // fPass==2 has the same triplet cuts as fPass==1
        }
        ++indx;
    }
    return fReaderTMVA2->EvaluateMVA("MLPBNN method");
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::InitTMVAtracks()
{
    // Initialize TMVA package for rejecting fake tracks (4- and 5-hits)
 
    // This loads the library
    // TMVA::Tools::Instance();
 
    // Default MVA methods to be trained + tested
    std::map<std::string, int> Use;
 
    Use["BDTD"] = 1;   // decorrelation + Adaptive Boost
    Float_t* var = fVarTMVA;
    fTmvaReaders.insert(pair<int, TMVA::Reader*>(-5, nullptr));
    fTmvaReaders.insert(pair<int, TMVA::Reader*>(-4, nullptr));
    fTmvaReaders.insert(pair<int, TMVA::Reader*>(4, nullptr));
    fTmvaReaders.insert(pair<int, TMVA::Reader*>(5, nullptr));
 
    // --- Create the Reader objects
 
    for (auto mit = fTmvaReaders.begin(); mit != fTmvaReaders.end(); ++mit) {
        TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
        int nhits = mit->first;
        fTmvaReaders[nhits] = reader;
 
        // Create a set of variables and declare them to the reader
        // The variable names MUST corresponds in name and type to those given in the weight file(s) used
 
        // FIXME
        int iii = 0;
        // reader->AddVariable( "sqrt(xb*xb+yb*yb)",                  &var[iii++] );
        reader->AddVariable("atan2(yb,xb)", &var[iii++]);
        // reader->AddVariable( "zb",                                 &var[iii++] );
        reader->AddVariable("tx", &var[iii++]);
        // reader->AddVariable( "ty",                                 &var[iii++] );
        reader->AddVariable("log(abs(momr))", &var[iii++]);
        reader->AddVariable("chi2", &var[iii++]);
        // reader->AddVariable( "log(sqrt(xextr*xextr+yextr*yextr))", &var[iii++] );
        // reader->AddVariable( "atan2(yextr,xextr)",                 &var[iii++] );
        reader->AddVariable("corrq[0]", &var[iii++]);
        reader->AddVariable("corrq[1]", &var[iii++]);
        reader->AddVariable("corrq[2]", &var[iii++]);
        reader->AddVariable("corrq[3]", &var[iii++]);
        if (TMath::Abs(nhits) == 5)
            reader->AddVariable("corrq[4]", &var[iii++]);
        reader->AddVariable("corrq3", &var[iii++]);
        // reader->AddVariable( "nmiss+rndm",                         &var[iii++] );
        reader->AddVariable("nseqmx", &var[iii++]);
        // reader->AddVariable( "leng",                               &var[iii++] );
        reader->AddVariable("corrn[0]", &var[iii++]);
        reader->AddVariable("corrn[1]", &var[iii++]);
        reader->AddVariable("corrn[2]", &var[iii++]);
        reader->AddVariable("corrn[3]", &var[iii++]);
        if (TMath::Abs(nhits) == 5)
            reader->AddVariable("corrn[4]", &var[iii++]);
        reader->AddVariable("nsigem", &var[iii++]);
 
        // --- Book the MVA methods
 
        TString dir = gSystem->Getenv("VMCWORKDIR");
        // TString dir    = "datasetTMVA1/weights";
        dir += "/parameters/reco/datasetTMVAV8/weights";
        // dir += nhits;
        dir += (nhits < 10) ? TMath::Abs(nhits) : nhits / 10;
        if (nhits < 0)
            dir += "neg";
        else
            dir += "pos";
        dir += "/";
        TString prefix = "TMVAClassification";
 
        // Book method(s)
        for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
            if (it->second) {
                TString methodName = TString(it->first) + TString(" method");
                TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
                reader->BookMVA(methodName, weightfile);
            }
        }
    }   // for (auto mit = fTmvaReaders.begin();
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::TMVAOutputTra()
{
    // Compute TMVA output for tracks
 
    // Fill cluster overlap map
    map<Int_t, set<Int_t>> mClusTr;
    int nTracks = fVectorArray->GetEntriesFast();
 
    for (int itr = 0; itr < nTracks; ++itr) {
        CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
        int nhit = tr->GetNStsHits();
        if (nhit < 4)
            continue;   // exclude 3-hit tracks
 
        for (int ih = 0; ih < nhit; ++ih) {
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(ih));
            // Check clusters
 
            for (Int_t side = 0; side < 2; ++side) {
                Int_t iclus = hit->GetDigi(side);
                if (mClusTr.count(iclus) == 0) {
                    set<Int_t> aaa;
                    mClusTr[iclus] = aaa;
                }
                mClusTr[iclus].insert(itr);
            }
        }
    }   // for (int itr = 0; itr < ntracks;
 
    Float_t* var = fVarTMVA;
 
    for (int itr = 0; itr < nTracks; ++itr) {
        CbmStsTrack* track = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
        int nhit = track->GetNStsHits();
        // if (nhit < 4 || nhit > 5) return 1.0;
        if (nhit < 4 || nhit > 5) {
            track->SetB(1.0);
            continue;
        }
        const FairTrackParam* param = track->GetParamFirst();
        Float_t momr = 1. / param->GetQp();
        int iq = (momr >= 0) ? 1 : -1;
        TVector3 posb, pose, mom3;
        Float_t corrq[19], corrn[19] = {0}, corrn3 = 0, corrq3 = 0;
        int nmiss = 0, ista0 = -1, nsigem = 0, nseqmx = 1, nseq = 1;
 
        for (int ih = 0; ih < nhit; ++ih) {
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(track->GetStsHitIndex(ih));
            Int_t ista = hit->GetStationNr() - 1;
            if (ih == 0) {
                hit->Position(posb);
                ista0 = ista;
            } else {
                if (ista - ista0 > 1) {
                    ++nmiss;
                    nseqmx = TMath::Max(nseqmx, nseq);
                    nseq = 1;
                } else
                    ++nseq;
                if (ih == nhit - 1)
                    hit->Position(pose);
            }
            corrq[ih] = hit->GetSignalDiv();
            corrq3 += TMath::Abs(corrq[ih]);
            ista0 = ista;
 
            for (Int_t side = 0; side < 2; ++side) {
                Int_t iclus = hit->GetDigi(side);
                int over = mClusTr[iclus].size();
                corrn[ih] += (over - 1);
            }
            corrn3 += corrn[ih];
            if (ista < 4)
                nsigem += 10;
            else
                nsigem += 1;
        }   // for (int ih = 0; ih < nhit;
 
        nseqmx = TMath::Max(nseqmx, nseq);
        pose -= posb;
        // Float_t leng = pose.Mag() / (nhit-1);
 
        // Extrapolate track to the target
        CbmKFTrack trKF = CbmKFTrack(*track);
        FairTrackParam par;
        // track.Extrapolate(xyzv[2]);
        trKF.Extrapolate(0.0);
        trKF.GetTrackParam(par);
        par.Momentum(mom3);
        momr = TMath::Sqrt(mom3.X() * mom3.X() + mom3.Z() * mom3.Z());
        par.Position(mom3);
 
        int iii = 0;
        // var[iii++] = TMath::Sqrt (xb*xb+yb*yb);
        // var[iii++] = posb.Pt();
        // var[iii++] = TMath::ATan2 (yb, xb);
        var[iii++] = posb.Phi();
        // var[iii++] = zb;
        // var[iii++] = posb.Z();
        var[iii++] = param->GetTx();
        var[iii++] = TMath::Log(TMath::Abs(momr));
        var[iii++] = track->GetChi2() / track->GetNDF();
        // var[iii++] = TMath::Log (TMath::Sqrt (xextr*xextr+yextr*yextr));
        // var[iii++] = TMath::ATan2 (yextr, xextr);
        // var[iii++] = TMath::Log (mom3.Pt());
        // var[iii++] = mom3.Phi();
        var[iii++] = corrq[0];
        var[iii++] = corrq[1];
        var[iii++] = corrq[2];
        var[iii++] = corrq[3];
        if (nhit == 5)
            var[iii++] = corrq[4];
        var[iii++] = corrq3;
        // var[iii++] = nmiss + gRandom->Rndm();
        var[iii++] = nseqmx;
        // var[iii++] = leng;
        // var[iii++] = corrn3;
        var[iii++] = corrn[0];
        var[iii++] = corrn[1];
        var[iii++] = corrn[2];
        var[iii++] = corrn[3];
        if (nhit == 5)
            var[iii++] = corrn[4];
        var[iii++] = nsigem;
 
        // return fTmvaReaders[nhit*iq]->EvaluateMVA("BDTD method");
        track->SetB(fTmvaReaders[nhit * iq]->EvaluateMVA("BDTD method"));
    }   // for (int itr = 0; itr < nTracks;
}
 
// -------------------------------------------------------------------------
 
void BmnStsVectorFinderV9::GoToTarget()
{
    // Evaluate vertex coordinates
 
    int ntracks = fVectorArray->GetEntriesFast();
    Double_t xv = 0, yv = 0, wwwxtot = 0, wwwytot = 0;
 
    for (int itr = 0; itr < ntracks; ++itr) {
        CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
        FairTrackParam param = *(tr->GetParamFirst());
        fitter.Extrapolate(&param, fXyzv[2], &param);
        // Double_t wwwx = 1 / TMath::Sqrt(param.GetCovariance(0,0));
        // Double_t wwwy = 1 / TMath::Sqrt(param.GetCovariance(1,1));
        Double_t wwwx = 1 / param.GetCovariance(0, 0);
        Double_t wwwy = 1 / param.GetCovariance(1, 1);
        xv += (param.GetX() * wwwx);
        yv += (param.GetY() * wwwy);
        wwwxtot += wwwx;
        wwwytot += wwwy;
        // cout << itr << " " << tr->GetNStsHits() << " " << param.GetX() << " " <<
        // TMath::Sqrt(param.GetCovariance(0,0)) << " "
        //    << param.GetY() << " " << TMath::Sqrt(param.GetCovariance(1,1)) << " " << tr->GetParamFirst()->GetZ() <<
        //    endl;
    }
    cout << " Vertex coords: " << fEvent << " " << fPass << " " << xv / wwwxtot << " " << yv / wwwytot << endl;
}
 
// -------------------------------------------------------------------------
// AZ-220626
 
void BmnStsVectorFinderV9::KillHitsTxP()
{
    // Kill hits on tracks in momentum space
 
    if (fHitConverter == nullptr) {
        fHitConverter = (BmnToCbmHitConverter*)FairRun::Instance()->GetTask("BMN to CBM Hits Converter");
        BmnStsVectorFinderV9* vf = (BmnStsVectorFinderV9*)FairRun::Instance()->GetTask("STS Vector Finder V9");
        // BmnStsVectorFinderV9 *vf = (BmnStsVectorFinderV9*) FairRun::Instance()->GetTask ("aaabbbccc");
        SetTracksToKill(vf->GetVectorArray());
    }
 
    int ntracks = fTracksToKill->GetEntriesFast(), nhits0 = fHitArray->GetEntriesFast(), nkilled = 0;
    FairTrackParam param;
    cout << " Initial ntracks: " << ntracks << endl;
 
    for (int itr = 0; itr < ntracks; ++itr) {
        CbmStsTrack* tr = (CbmStsTrack*)fTracksToKill->UncheckedAt(itr);
        CbmKFTrack kftr = CbmKFTrack(*tr, kTRUE);
        kftr.Extrapolate(0.0);   // extrapolate to the target
        kftr.GetTrackParam(param);
        Double_t tx = param.GetTx();
        Double_t ty = param.GetTy();
        Double_t pmom = 1 / param.GetQp();
        int iq = (pmom >= 0) ? 1 : -1;
        int nhits = tr->GetNStsHits();
 
        for (int ih = 0; ih < nhits; ++ih) {
            int ihh = tr->GetStsHitIndex(ih);
            CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ihh);
            if (hit == nullptr)
                continue;
            Bool_t ok = fHitConverter->CorrectEffic(hit, iq, tx, ty, TMath::Abs(pmom), 1);
            // if (!ok) fHitArray->RemoveAt (ihh);
            if (!ok) {
                hit->SetUniqueID(999);   // exclude hit from tracking
                ++nkilled;
            }
        }
    }
    fHitArray->Compress();
    cout << " !!! Hits: " << nhits0 << " " << fHitArray->GetEntriesFast() << " " << nkilled << endl;
}
 
// -------------------------------------------------------------------------
 
ClassImp(BmnStsVectorFinderV9)