BmnVertexFinder.cxx

Go to the documentation of this file.

#include "BmnVertexFinder.h"
 
#include "BmnMath.h"
#include "TCanvas.h"
#include "TFitResult.h"
#include "TGraph.h"
#include "vector"
 
using namespace std;
using namespace TMath;
 
BmnVertexFinder::BmnVertexFinder(Int_t period) {
    fPeriodId = period;
    fEventNo = 0;
    fGlobalTracksArray = NULL;
    fRobustRefit = kTRUE;
    fGlobalTracksBranchName = "BmnGlobalTrack";
    fVertexBranchName = "BmnVertex";
}
 
BmnVertexFinder::~BmnVertexFinder() {
    delete fKalman;
}
 
InitStatus BmnVertexFinder::Init() {
    if (fVerbose > 1) cout << "=========================== Vertex finder init started ====================" << endl;
 
    fKalman = new BmnKalmanFilter();
 
    //Get ROOT Manager
    FairRootManager* ioman = FairRootManager::Instance();
    if (NULL == ioman) Fatal("Init", "FairRootManager is not instantiated");
 
    fGlobalTracksArray = (TClonesArray*)ioman->GetObject(fGlobalTracksBranchName);  //in
    if (!fGlobalTracksArray) {
        cout << "BmnVertexFinder::Init(): branch " << fGlobalTracksBranchName << " not found! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
 
    fVertexArray = new TClonesArray("BmnVertex", 1);  //out
    ioman->Register(fVertexBranchName, "GEM", fVertexArray, kTRUE);
 
    if (fVerbose > 1) cout << "=========================== Vertex finder init finished ===================" << endl;
 
    return kSUCCESS;
}
 
void BmnVertexFinder::Exec(Option_t* opt) {
    TStopwatch sw;
    sw.Start();
 
    if (!IsActive())
        return;
 
    if (fVerbose > 1) cout << "======================== Vertex finder exec started  ======================" << endl;
    if (fVerbose > 1) cout << "Event number: " << fEventNo++ << endl;
 
    fVertexArray->Delete();
 
    new ((*fVertexArray)[fVertexArray->GetEntriesFast()]) BmnVertex(FindPrimaryVertex(fGlobalTracksArray));
    new ((*fVertexArray)[fVertexArray->GetEntriesFast()]) BmnVertex(FindSecondaryVertex(fGlobalTracksArray));
 
    if (fVerbose > 1) cout << "\n======================== Vertex finder exec finished ======================" << endl;
 
    sw.Stop();
    fTime += sw.RealTime();
}
 
BmnVertex BmnVertexFinder::FindPrimaryVertex(TClonesArray* tracks) {
    //cout << "Primary" << endl;
    const Double_t kRange = (fPeriodId == 8) ? 200.0 : 50.0;
    const Double_t xVertMin = (fPeriodId == 8) ? -10.0 : -0.5;
    const Double_t xVertMax = (fPeriodId == 8) ? +10.0 : +2.0;
    const Double_t yVertMin = (fPeriodId == 8) ? -10.0 : -5.0;
    const Double_t yVertMax = (fPeriodId == 8) ? +10.0 : -2.5;
    Double_t roughVertexZ = (fPeriodId == 8) ? 0.0 : -1.0;
    Int_t nPrim = 0;
 
    for (Int_t iTr = 0; iTr < tracks->GetEntriesFast(); ++iTr) {
        BmnGlobalTrack* track = (BmnGlobalTrack*)tracks->At(iTr);
        FairTrackParam par0 = *(track->GetParamFirst());
        if (fKalman->TGeoTrackPropagate(&par0, roughVertexZ, (par0.GetQp() > 0.) ? 2212 : -211, NULL, NULL) == kBMNERROR) {
            continue;
        }
 
        if (par0.GetX() > xVertMin && par0.GetX() < xVertMax && par0.GetY() > yVertMin && par0.GetY() < yVertMax) {
            track->SetFlag(0);  //primary track
            nPrim++;
        } else {
            track->SetFlag(1);  //secondary track
        }
    }
 
    if (nPrim < 2) {
        return BmnVertex();
    }
 
    Double_t vz = FindVZByVirtualPlanes(roughVertexZ, kRange, tracks, 0);
    if (Abs(vz - roughVertexZ) > kRange) {
        for (Int_t iTr = 0; iTr < tracks->GetEntriesFast(); ++iTr) {
            BmnGlobalTrack* track = (BmnGlobalTrack*)tracks->At(iTr);
            track->SetFlag(1);
        }
        return BmnVertex();
    }
 
    vector<Double_t> xHits;
    vector<Double_t> yHits;
    vector<Int_t> indexes;
    for (Int_t iTr = 0; iTr < tracks->GetEntriesFast(); ++iTr) {
        BmnGlobalTrack* track = (BmnGlobalTrack*)tracks->At(iTr);
        if (track->GetFlag() != 0) continue;
        FairTrackParam par0 = *(track->GetParamFirst());
        if (fKalman->TGeoTrackPropagate(&par0, vz, (par0.GetQp() > 0.) ? 2212 : -211, NULL, NULL) == kBMNERROR) {
            continue;
        }
        xHits.push_back(par0.GetX());
        yHits.push_back(par0.GetY());
        indexes.push_back(iTr);
    }
 
    Double_t vx = Mean(xHits.begin(), xHits.end());
    Double_t vy = Mean(yHits.begin(), yHits.end());
 
    Double_t rRMS = CalcRms2D(xHits, yHits);
 
    return BmnVertex(vx, vy, vz, rRMS, 0, xHits.size(), TMatrixFSym(3), -1, vector<Int_t>());
}
 
BmnVertex BmnVertexFinder::FindSecondaryVertex(TClonesArray* tracks) {
    //cout << "Sacondary" << endl;
    Double_t roughVertexZ = (fPeriodId == 8) ? 0.0 : (fPeriodId == 7) ? -1.0 : (fPeriodId == 6) ? -21.9 : 0.0;
    const Double_t kRange = 300.0;
 
    Int_t nSec = 0;
 
    for (Int_t iTr = 0; iTr < tracks->GetEntriesFast(); ++iTr) {
        BmnGlobalTrack* track = (BmnGlobalTrack*)tracks->At(iTr);
        if (track->GetFlag() == 1) nSec++;  //secondary track
    }
 
    if (nSec < 2) {
        return BmnVertex();
    }
 
    Double_t vz = FindVZByVirtualPlanes(roughVertexZ, kRange, tracks, 1);
 
    vector<Double_t> xHits;
    vector<Double_t> yHits;
    vector<Int_t> indexes;
    for (Int_t iTr = 0; iTr < tracks->GetEntriesFast(); ++iTr) {
        BmnGlobalTrack* track = (BmnGlobalTrack*)tracks->At(iTr);
        if (track->GetFlag() != 1) continue;
        FairTrackParam par0 = *(track->GetParamFirst());
        if (fKalman->TGeoTrackPropagate(&par0, vz, (par0.GetQp() > 0.) ? 2212 : -211, NULL, NULL) == kBMNERROR) {
            continue;
        }
        xHits.push_back(par0.GetX());
        yHits.push_back(par0.GetY());
        indexes.push_back(iTr);
    }
 
    Double_t vx = Mean(xHits.begin(), xHits.end());
    Double_t vy = Mean(yHits.begin(), yHits.end());
 
    Double_t rRMS = CalcRms2D(xHits, yHits);
 
    return BmnVertex(vx, vy, vz, rRMS, 0, xHits.size(), TMatrixFSym(3), -1, vector<Int_t>());
}
 
Float_t BmnVertexFinder::FindVZByVirtualPlanes(Float_t z_0, Float_t range, TClonesArray* tracks, Int_t flag) {
    const Int_t nPlanes = 3;
    Float_t minZ = z_0;
 
    while (range >= 0.01) {
        Float_t zMax = minZ + range;
        Float_t zMin = minZ - range;
        Float_t zStep = (zMax - zMin) / (nPlanes - 1);
 
        vector<Double_t> xHits[nPlanes];
        vector<Double_t> yHits[nPlanes];
        Float_t zPlane[nPlanes];
        Float_t rRMS[nPlanes] = {0};
 
        for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane)
            zPlane[iPlane] = zMax - iPlane * zStep;
 
        Int_t nOkTr = 0;
 
        for (Int_t iTr = 0; iTr < tracks->GetEntriesFast(); ++iTr) {
            BmnGlobalTrack* track = (BmnGlobalTrack*)tracks->At(iTr);
            if (track->GetFlag() != flag) continue;
            FairTrackParam par0 = *(track->GetParamFirst());
            //cout << "tr" << iTr << endl;
            Double_t xTr[nPlanes];
            Double_t yTr[nPlanes];
            Bool_t trOk = kTRUE;
            for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
                if (fKalman->TGeoTrackPropagate(&par0, zPlane[iPlane], (par0.GetQp() > 0.) ? 2212 : -211, NULL, NULL) == kBMNERROR) {
                    track->SetFlag(-1);  //maybe 666???
                    trOk = kFALSE;
                    break;
                }
 
                //cout << " " << par0.GetX() << " " << par0.GetY() << " " << par0.GetZ() << endl;
                xTr[iPlane] = par0.GetX();
                yTr[iPlane] = par0.GetY();
                //xHits[iPlane].push_back(par0.GetX());
                //yHits[iPlane].push_back(par0.GetY());
            }
 
            if (trOk) {
                nOkTr++;
                for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
                    xHits[iPlane].push_back(xTr[iPlane]);
                    yHits[iPlane].push_back(yTr[iPlane]);
                }
            }
        }
 
        if (nOkTr < 2) return -1000.0;
 
        //Double_t minRMS = DBL_MAX;
 
        //Calculation minZ as minimal value
        // for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
        //     rRMS[iPlane] = CalcMeanDist(xHits[iPlane], yHits[iPlane]);
        //     if (rRMS[iPlane] < minRMS) {
        //         minRMS = rRMS[iPlane];
        //         minZ = zPlane[iPlane];
        //     }
        // }
 
        //Calculation minZ as minimum of parabola
        for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
            rRMS[iPlane] = CalcMeanDist(xHits[iPlane], yHits[iPlane]);
        }
        TGraph* vertex = new TGraph(nPlanes, zPlane, rRMS);
        TFitResultPtr ptr = vertex->Fit("pol2", "QFS");
        Float_t b = ptr->Parameter(1);
        Float_t a = ptr->Parameter(2);
        minZ = -b / (2 * a);
        delete vertex;
 
        range /= 2;
 
        // for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane)
        //     //rRMS[iPlane] = CalcRms2D(xHits[iPlane], yHits[iPlane]);
        // rRMS[iPlane] = CalcMeanDist(xHits[iPlane], yHits[iPlane]);
 
        //TCanvas* c = new TCanvas("c1", "c1", 1200, 800);
        // //c->Divide(1, 1);
        // //c->cd(1);
        //TGraph* vertex = new TGraph(nPlanes, zPlane, rRMS);
        // TFitResultPtr ptr = vertex->Fit("pol2", "QFS");
        // //        TF1 *fit_func = vertex->GetFunction("pol2");
        // Float_t b = ptr->Parameter(1);
        // Float_t a = ptr->Parameter(2);
        // minZ = -b / (2 * a);
        // range /= 2;
        //vertex->Draw("AP*");
        //c->SaveAs("tmp.pdf");
        //getchar();
        //delete vertex;
    }
    return minZ;
}
 
void BmnVertexFinder::Finish() {
    ofstream outFile;
    outFile.open("QA/timing.txt", ofstream::app);
    outFile << "Vertex Finder Time: " << fTime;
    if (fVerbose == 0) cout << "Work time of the GEM vertex finder: " << fTime << endl;
}
 
Double_t BmnVertexFinder::CalcRms2D(vector<Double_t> x, vector<Double_t> y) {
    Double_t xMean = 0.0;
    Double_t yMean = 0.0;
    Double_t rms = 0.0;
    for (size_t iHit = 0; iHit < x.size(); ++iHit) {
        xMean += x.at(iHit);
        yMean += y.at(iHit);
    }
    xMean /= x.size();
    yMean /= y.size();
 
    for (size_t iHit = 0; iHit < x.size(); ++iHit) {
        rms += (Sq(x.at(iHit) - xMean) + Sq(y.at(iHit) - yMean));
    }
 
    rms /= x.size();
    return Sqrt(rms);
}
 
Double_t BmnVertexFinder::CalcMeanDist(vector<Double_t> x, vector<Double_t> y) {
    Double_t sumDist = 0.0;
    Int_t nPairs = 0;
    for (size_t i = 0; i < x.size(); ++i) {
        for (size_t j = i + 1; j < x.size(); ++j) {
            sumDist += Sqrt(Sq(x[i] - x[j]) + Sq(y[i] - y[j]));
            nPairs++;
        }
    }
    return sumDist / nPairs;  // calc. ave. dist value
}