BmnEfficiency.h

Go to the documentation of this file.

#ifndef BMNEFF_H
#define BMNEFF_H 1
 
#include "BmnNewFieldMap.h"
 
#include <BmnGlobalTrack.h>
#include <BmnHit.h>
#include <BmnInnerTrackerGeometryDraw.h>
#include <BmnKalmanFilter.h>
#include <BmnVertex.h>
#include <CbmVertex.h>
#include <DstEventHeader.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TEfficiency.h>
#include <TGeoManager.h>
#include <TH1.h>
#include <TH2.h>
#include <TLine.h>
#include <TNamed.h>
#include <TObject.h>
#include <TStyle.h>
#include <UniRun.h>
#include <vector>
 
using namespace std;
 
// A class to be used as a store when calculating efficiency ...
 
class EffStore : public TObject
{
  public:
    EffStore()
        : detector("")
        , station(-1)
        , module(-1)
        , zone("")
        , nominator(0)
        , denominator(0)
    {}
 
    EffStore(TString det, Int_t stat, Int_t mod, TString z = "")
        : detector(det)
        , station(stat)
        , module(mod)
        , zone(z)
        , nominator(0)
        , denominator(0)
    {
        hitPositionsAndEff.resize(0);
    }
 
    EffStore(TString det, Int_t stat, TString z = "")
        : detector(det)
        , station(stat)
        , module(-1)
        , zone(z)
        , nominator(0)
        , denominator(0)
    {
        hitPositionsAndEff.resize(0);
    }
 
    void IncreaseNominatorByUnity() { nominator++; }
 
    void IncreaseDenominatorByUnity() { denominator++; }
 
    void AddHitCoordinates(BmnHit* hit, Double_t eff = 1.)
    {
        hitPositionsAndEff.push_back(TLorentzVector(hit->GetX(), hit->GetY(), hit->GetZ(), eff));
    }
 
    Double_t Efficiency()
    {
        Double_t eff = 0.;
        if (denominator != 0)
            eff = 1. * nominator / denominator;
 
        return eff;
    }
 
    TString Detector() { return detector; }
    Int_t Station() { return station; }
    Int_t Module() { return module; }
    TString Zone() { return zone; }
 
    vector<TLorentzVector> CoordinatesAndEffs() { return hitPositionsAndEff; }
 
    virtual ~EffStore() {}
 
  private:
    TString detector;   // GEM or SILICON
    Int_t station;
    Int_t module;
    TString zone;   // big or hot (valid for GEM only!)
 
    Int_t nominator;
    Int_t denominator;
 
    vector<TLorentzVector> hitPositionsAndEff;
 
    ClassDef(EffStore, 1)
};
 
class inputForEfficiency
{
  public:
    inputForEfficiency()
        : det("GEM")
        , nGemCut(4)
        , nSilCut(2)
    {
        eContainers = nullptr;
    }
 
    virtual ~inputForEfficiency() {}
 
    // Move assignment
 
    inputForEfficiency& operator=(inputForEfficiency&& input)
    {
        if (&input == this)
            return *this;
 
        det = input.det;
        nGemCut = input.nGemCut;
        nSilCut = input.nSilCut;
 
        sfIndices.first = input.sfIndices.first;
        sfIndices.second = input.sfIndices.second;
 
        flParams.first = input.flParams.first;
        flParams.second = input.flParams.second;
 
        gPairMap = input.gPairMap;
        sPairMap = input.sPairMap;
 
        eContainers = input.eContainers;
 
        for (auto it : input.usedStats)
            usedStats.push_back(it);
 
        for (auto it : input.hits)
            hits.push_back(it);
 
        input.hits.clear();
        input.usedStats.clear();
 
        if (input.eContainers)
            input.eContainers = nullptr;
 
        return *this;
    }
 
    // Detector specified, hit arrays and stats used ...
    TString det;
    vector<BmnHit*> hits;
    vector<Bool_t> usedStats;
 
    // Indices and track params ....
    pair<Int_t, Int_t> sfIndices;
    pair<FairTrackParam, FairTrackParam> flParams;
 
    TClonesArray* eContainers;
 
    // Y-ranges ...
    map<Int_t, vector<pair<Double_t, Double_t>>> gPairMap;   // stat --> corresponding set of Y-pairs
    map<Int_t, vector<pair<Double_t, Double_t>>> sPairMap;
 
    // hit cuts ...
    Int_t nGemCut;
    Int_t nSilCut;
};
 
class InnerTrackerParams : public TObject
{
  public:
    InnerTrackerParams() {}
 
    InnerTrackerParams(Int_t stat, Double_t xMin, Double_t xMax, Double_t yMin, Double_t yMax)
        : fStation(stat)
        , fModule(-1)
        , fXMin(xMin)
        , fXMax(xMax)
        , fYMin(yMin)
        , fYMax(yMax)
    {}
 
    InnerTrackerParams(Int_t stat, Int_t mod, Double_t xMin, Double_t xMax, Double_t yMin, Double_t yMax)
        : fStation(stat)
        , fModule(mod)
        , fXMin(xMin)
        , fXMax(xMax)
        , fYMin(yMin)
        , fYMax(yMax)
    {}
 
    InnerTrackerParams(Int_t stat, Int_t mod, Double_t x, Double_t y, Double_t z)
        : fStation(stat)
        , fModule(mod)
    {
        hitX = x;
        hitY = y;
        hitZ = z;
    }
 
    virtual ~InnerTrackerParams() { ; }
 
    Int_t station() { return fStation; }
 
    Int_t module() { return fModule; }
 
    pair<Double_t, Double_t> x() { return make_pair(fXMin, fXMax); }
 
    pair<Double_t, Double_t> y() { return make_pair(fYMin, fYMax); }
 
    void SetZ(Double_t z) { fZ = z; }
 
    Double_t z() { return fZ; }
 
    TVector3 xyz() { return TVector3(hitX, hitY, hitZ); }
 
  protected:
    Int_t fStation;
    Int_t fModule;
 
    Double_t fXMin;
    Double_t fXMax;
 
    Double_t fYMin;
    Double_t fYMax;
 
    Double_t fZ;   // position of station (at given module if defined) ...
 
    //
    Double_t hitX;
    Double_t hitY;
    Double_t hitZ;
 
    ClassDef(InnerTrackerParams, 1)
};
 
class EffStore2D : public EffStore
{
  public:
    EffStore2D() { ; }
 
    EffStore2D(TString det, Int_t stat, map<Int_t, vector<pair<Double_t, Double_t>>> pairs)
        : EffStore(det, stat)
        , hEffGem(nullptr)
        , hEffGem1d(nullptr)
        , hEffSil(nullptr)
        , hEffSil1d(nullptr)
    {
 
        const Int_t nQp = 2;   // positive and negative ...
 
        enum QpValue
        {
            pos,
            neg
        };
 
        hEffGem = new TEfficiency*[nQp];
        hEffSil = new TEfficiency*[nQp];
        hEffSilXP = new TEfficiency*[nQp];
 
        if (det.Contains("GEM")) {
            Int_t nBinsX = 160.;   // 80.; // 160;
            Int_t nBinsY = 45;
 
            pair<Double_t, Double_t> y = make_pair(0., 45.);
 
            Double_t xLow = -80.;
            Double_t xUp = +80.;
 
            for (Int_t iQp = 0; iQp < nQp; iQp++) {
                TString Qp = (iQp == QpValue::pos) ? "Q_{p} > 0" : "Q_{p} < 0";
 
                hEffGem[iQp] = new TEfficiency(Form("Detector# %s, Station# %d, %s", det.Data(), stat, Qp.Data()),
                                               Form("Stat %d, %s; X [cm]; Y [cm]", stat, Qp.Data()), nBinsX, xLow, xUp,
                                               nBinsY, y.first, y.second);
            }
 
            // Preparing 1D-efficiencies (x for different y-ranges) ...
            if (pairs.size()) {
                nBinsX = 20.;   // 80.; // 160;
 
                hEffGem1d = new TEfficiency**[nQp];
 
                const Int_t nYRanges = pairs.find(stat)->second.size();
 
                for (Int_t iQp = 0; iQp < nQp; iQp++) {
                    hEffGem1d[iQp] = new TEfficiency*[nYRanges];
 
                    TString Qp = (iQp == QpValue::pos) ? "Q_{p} > 0" : "Q_{p} < 0";
 
                    for (Int_t iRange = 0; iRange < nYRanges; iRange++) {
                        Double_t y0 = pairs.find(stat)->second.at(iRange).first;
                        Double_t y1 = pairs.find(stat)->second.at(iRange).second;
 
                        hEffGem1d[iQp][iRange] = new TEfficiency(
                            Form("Detector# %s, Station# %d, (%G < y < %G [cm]), %s", det.Data(), stat, y0, y1,
                                 Qp.Data()),
                            Form("Stat %d, (%G < y < %G [cm]), %s; X [cm]; #varepsilon", stat, y0, y1, Qp.Data()),
                            nBinsX, xLow, xUp);
                    }
                }
            }
        } else {
            Int_t nBinsX = 80;
            Int_t nBinsY = 45;
 
            TClonesArray* sil = new TClonesArray("InnerTrackerParams");
 
            // Adding necessary information ...
            // stat 0
            new ((*sil)[sil->GetEntriesFast()]) InnerTrackerParams(0, -6., +6., -5., +8.);
 
            // stat 1
            new ((*sil)[sil->GetEntriesFast()]) InnerTrackerParams(1, -6., +6., -5., +8.);
 
            // stat 2
            new ((*sil)[sil->GetEntriesFast()]) InnerTrackerParams(2, -12., +12., -11., +14.);
 
            for (Int_t iQp = 0; iQp < nQp; iQp++) {
                TString Qp = (iQp == QpValue::pos) ? "Q_{p} > 0" : "Q_{p} < 0";
                hEffSil[iQp] = new TEfficiency(Form("Detector# %s, Station# %d, %s", det.Data(), stat, Qp.Data()),
                                               Form("Stat %d, %s; X [cm]; Y [cm]", stat, Qp.Data()), nBinsX,
                                               getRanges(sil, stat)[0], getRanges(sil, stat)[1], nBinsY,
                                               getRanges(sil, stat)[2], getRanges(sil, stat)[3]);
 
                nBinsY = 55;
 
                hEffSilXP[iQp] =
                    new TEfficiency(Form("Detector# %s, Station# %d (xP), %s", det.Data(), stat, Qp.Data()),
                                    Form("Stat %d, %s; X [cm]; P [GeV/c]", stat, Qp.Data()), nBinsX,
                                    getRanges(sil, stat)[0], getRanges(sil, stat)[1], nBinsY, 0., 5.5);
            }
 
            // Preparing 1D-efficiencies (x for different y-ranges) ...
            if (pairs.size()) {
                nBinsX = 20.;   // 80.; // 160;
 
                hEffSil1d = new TEfficiency**[nQp];
 
                const Int_t nYRanges = pairs.find(stat)->second.size();
 
                for (Int_t iQp = 0; iQp < nQp; iQp++) {
                    hEffSil1d[iQp] = new TEfficiency*[nYRanges];
 
                    TString Qp = (iQp == QpValue::pos) ? "Q_{p} > 0" : "Q_{p} < 0";
 
                    for (Int_t iRange = 0; iRange < nYRanges; iRange++) {
                        Double_t y0 = pairs.find(stat)->second.at(iRange).first;
                        Double_t y1 = pairs.find(stat)->second.at(iRange).second;
 
                        hEffSil1d[iQp][iRange] = new TEfficiency(
                            Form("Detector# %s, Station# %d, (%G < y < %G [cm]), %s", det.Data(), stat, y0, y1,
                                 Qp.Data()),
                            Form("Stat %d, (%G < y < %G [cm]), %s; X [cm]; #varepsilon", stat, y0, y1, Qp.Data()),
                            nBinsX, getRanges(sil, stat)[0], getRanges(sil, stat)[1]);
                    }
                }
            }
 
            delete sil;
        }
    }
 
    TEfficiency* efficiency(Int_t qpBin = 0, TString hType = "", TString det = "GEM")
    {
        // qpBin = 0 (positive), 1 (negative)
        if (hType.IsNull())
            return det.Contains("GEM") ? hEffGem[qpBin] : hEffSil[qpBin];
        else if (det.Contains("SILICON") && !hType.IsNull())
            return hEffSilXP[qpBin];
        else
            return nullptr;
    }
 
    TEfficiency* efficiency4range(Int_t range, Int_t qpBin = 0, TString det = "GEM")
    {
        // qpBin = 0 (positive), 1 (negative)
 
        if (det.Contains("GEM") && hEffGem1d[qpBin][range])
            return hEffGem1d[qpBin][range];
        else if (det.Contains("SIL") && hEffSil1d[qpBin][range])
            return hEffSil1d[qpBin][range];
 
        else
            return nullptr;
    }
 
    virtual ~EffStore2D()
    {
        // Should be implemented further ...
    }
 
  private:
    TEfficiency** hEffGem;
    TEfficiency*** hEffGem1d;
 
    TEfficiency** hEffSil;      // epslion(x, y)
    TEfficiency** hEffSilXP;    // epsilon(x, P)
    TEfficiency*** hEffSil1d;   // epsilon(x) at different y-slices
 
  private:
    vector<Double_t> getRanges(TClonesArray* sil, Int_t stat)
    {
 
        vector<Double_t> tmp;
 
        for (Int_t iEntry = 0; iEntry < sil->GetEntriesFast(); iEntry++) {
            InnerTrackerParams* silEle = (InnerTrackerParams*)sil->UncheckedAt(iEntry);
 
            if (silEle->station() == stat) {
                tmp.push_back(silEle->x().first);
                tmp.push_back(silEle->x().second);
                tmp.push_back(silEle->y().first);
                tmp.push_back(silEle->y().second);
 
                break;
            }
        }
 
        return tmp;
    }
 
    ClassDef(EffStore2D, 1)
};
 
// A class to be used as a store for calculated residuals ...
 
class Residuals : public EffStore
{
  public:
    Residuals()
        : detector("")
        , station(-1)
        , module(-1)
    {
        ;
    }
 
    Residuals(TString det, Int_t stat, Int_t mod)
        : EffStore(det, stat, mod)
        , detector(det)
        , station(stat)
        , module(mod)
    {
        ;
    }
 
    virtual ~Residuals() { ; }
 
    void SetXY(Double_t xR, Double_t yR)
    {
        x.push_back(xR);
        y.push_back(yR);
    }
 
    vector<Double_t> GetX() { return x; }
 
    vector<Double_t> GetY() { return y; }
 
  private:
    TString detector;   // GEM or SILICON
    Int_t station;
    Int_t module;
 
    vector<Double_t> x;
    vector<Double_t> y;
 
    ClassDef(Residuals, 1)
};
 
class BmnEfficiency : public TNamed
{
  public:
    BmnEfficiency();
    BmnEfficiency(FairRunAna*, BmnInnerTrackerGeometryDraw*, TString, Int_t nEvents = 0);
    BmnEfficiency(FairRunAna*, TString, Int_t nEvents = 0);
 
    void SetMinNHitsPerGlobTrack(Int_t nHits) { fNHits = nHits; }
 
    void SetMinNHitsPerSiliconTrack(Int_t nHits) { fNHitsSilicon = nHits; }
 
    void SetMinNHitsPerGemTrack(Int_t nHits) { fNHitsGem = nHits; }
 
    void Efficiency(TClonesArray*,
                    TClonesArray* effSilicon,
                    map<Int_t, vector<pair<Double_t, Double_t>>>,
                    map<Int_t, vector<pair<Double_t, Double_t>>>);
 
    virtual ~BmnEfficiency()
    {
        if (fMagField)
            delete fMagField;
        if (dstChain)
            delete dstChain;
        if (fKalman)
            delete fKalman;
    }
 
    TString GetGemZone(BmnHit*);
 
    // Overloading =  ...
 
    BmnEfficiency& operator=(const BmnEfficiency& eff)
    {
        // Check assignement to itself
        if (this == &eff)
            return *this;
 
        dstChain = eff.dstChain;
        fVertices = eff.fVertices;
        fGlobTracks = eff.fGlobTracks;
 
        fSiliconTracks = eff.fSiliconTracks;
        fSiliconHits = eff.fSiliconHits;
 
        fGemTracks = eff.fGemTracks;
        fGemHits = eff.fGemHits;
 
        fMagField = eff.fMagField;
        fField = eff.fField;
        fKalman = eff.fKalman;
 
        return *this;
    }
 
  private:
    Bool_t isGoodDst;
    // Acceptance maps to be used ...
    map<pair<Int_t, Int_t>, Double_t> fStatZ;                    // (stat, mod) --> Z
    map<pair<Int_t, Int_t>, vector<Double_t>> fStatAcceptance;   // (stat, mod) --> vector of boarders ...
    Int_t fNHitsSilicon;
    Int_t fNHitsGem;
 
    void fillEfficiency(inputForEfficiency&& in);
    void FillAcceptanceMaps(TClonesArray*, TClonesArray* SILs = nullptr);
 
    inline TString GetDetector(BmnHit* hit)
    {
        Double_t z = hit->GetZ();
 
        TString det = "";
        if (z > gem->GetStation(0)->GetZPosition() - 5.)
            det = "GEM";
        else
            det = "SILICON";
 
        return det;
    }
 
    inline Double_t FindZ(Int_t stat) { return fStatZ.find(make_pair(stat, 0))->second; }
 
    Bool_t isVirtualHitInAcceptance(Int_t,
                                    vector<BmnHit*>,
                                    pair<FairTrackParam, FairTrackParam>,
                                    pair<Double_t, Double_t>&);
 
    Int_t FindBin(Double_t y, Int_t stat, map<Int_t, vector<pair<Double_t, Double_t>>> pairMap)
    {
        Int_t bin = -1;
 
        for (size_t iRange = 0; iRange < pairMap.find(stat)->second.size(); iRange++) {
            Double_t min = pairMap.find(stat)->second.at(iRange).first;
            Double_t max = pairMap.find(stat)->second.at(iRange).second;
 
            if (y > min && y < max) {
                bin = iRange;
                break;
            }
        }
 
        return bin;
    }
 
  protected:
    // Geometries
    BmnGemStripStationSet* gem;
    BmnSiliconStationSet* silicon;
 
    DstEventHeader* fHeader;
 
    TClonesArray* fInnerHits;
    TClonesArray* fGemHits;
    TClonesArray* fSiliconHits;
    TClonesArray* fGlobTracks;
    TClonesArray* fGemTracks;
    TClonesArray* fSiliconTracks;
    TClonesArray* fVertices;
 
    TChain* dstChain;
    Int_t fNEvents;
 
    FairField* fField;
    BmnNewFieldMap* fMagField;
    BmnKalmanFilter* fKalman;
 
    Int_t fNHits;
 
    ClassDef(BmnEfficiency, 1)
};
 
#endif