BmnDedx.cxx

Go to the documentation of this file.

#include "BmnDedx.h"
 
BmnDedx::BmnDedx()
    : fMinMom(0.0)
    , fMaxMom(20.0)
    , fNormDedxThresh(2.8)
    , fNSt(11)
    , fNSect({6, 10, 14, 18, 8, 8, 8, 8, 8, 8, 8})
    , fNMomSlices(6)
    , fMinMomBound({0.0, 0.5, 1.0, 1.5, 2.0, 3.0})
    , fMaxMomBound({0.5, 1.0, 2.0, 3.0, 5.0, 20.0})
    , fRunId(-1)
    , fStsTracks(nullptr)
    , fStsHits(nullptr)
    , fLowGemClusters(nullptr)
    , fLowFsdClusters(nullptr)
    , fGlobTracks(nullptr)
    , fDedx(0.0)
    , fTime(0.0)
{
    const Float_t fpar[5] = {-295.415, -2.71164, -0.0606733, 1.93517, 3.99467e-6};
    fAlephFit_He4 =
        new TF1("alephFit_He4", "[0]/(x^([3]))*([1]-x^([3])-TMath::Log([2]+1/(x^([4]))))", fMinMom, fMaxMom);
    fAlephFit_He4->SetParameters(fpar[0], fpar[1], fpar[2], fpar[3], fpar[4]);
}
 
BmnDedx::BmnDedx(Int_t runid)
    : fMinMom(0.0)
    , fMaxMom(20.0)
    , fNormDedxThresh(2.8)
    , fNSt(11)
    , fNSect({6, 10, 14, 18, 8, 8, 8, 8, 8, 8, 8})
    , fNMomSlices(6)
    , fMinMomBound({0.0, 0.5, 1.0, 1.5, 2.0, 3.0})
    , fMaxMomBound({0.5, 1.0, 2.0, 3.0, 5.0, 20.0})
    , fRunId(runid)
    , fStsTracks(nullptr)
    , fStsHits(nullptr)
    , fLowGemClusters(nullptr)
    , fLowFsdClusters(nullptr)
    , fGlobTracks(nullptr)
    , fDedx(0.0)
    , fTime(0.0)
{
    const Float_t fpar[5] = {-295.415, -2.71164, -0.0606733, 1.93517, 3.99467e-6};
    fAlephFit_He4 =
        new TF1("alephFit_He4", "[0]/(x^([3]))*([1]-x^([3])-TMath::Log([2]+1/(x^([4]))))", fMinMom, fMaxMom);
    fAlephFit_He4->SetParameters(fpar[0], fpar[1], fpar[2], fpar[3], fpar[4]);
}
 
void BmnDedx::Exec(Option_t* opt)
{
 
    TStopwatch sw;
    sw.Start();
 
    if (!IsActive())
        return;
 
    if (fVerbose > 1)
        cout << "\n======================== BmnDedx::Exec started =====================\n" << endl;
 
    for (Int_t iTr = 0; iTr < fGlobTracks->GetEntriesFast(); ++iTr) {
        BmnGlobalTrack* tr = (BmnGlobalTrack*)fGlobTracks->At(iTr);
        CalcDedx(tr);
    }
 
    sw.Stop();
    fTime += sw.RealTime();
 
    if (fVerbose > 1)
        cout << "\n======================== BmnDedx::Exec started ====================\n" << endl;
}
 
void BmnDedx::Finish()
{
    printf("Work time of BmnDedx: %4.2f sec.\n", fTime);
}
 
InitStatus BmnDedx::Init()
{
    if (fVerbose > 1)
        cout << "\n======================== BmnDedx::Init finished ====================\n" << endl;
 
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman)
        Fatal("Init", "FairRootManager is not instantiated");
 
    fStsTracks = (TClonesArray*)ioman->GetObject("StsVector");
    if (!fStsTracks) {
        cout << "BmnDedx::Init(): branch StsVector not found! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
    fStsHits = (TClonesArray*)ioman->GetObject("StsHit");
    if (!fStsHits) {
        cout << "BmnDedx::Init(): branch StsHit not found! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
    fGlobTracks = (TClonesArray*)ioman->GetObject("BmnGlobalTrack");
    if (!fGlobTracks) {
        cout << "BmnDedx::Init(): branch BmnGlobalTrack not found! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
    fLowGemClusters = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
    if (!fLowGemClusters) {
        cout << "BmnDedx::Init(): branch BmnGemLowerCluster not found! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
    fLowFsdClusters = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
    if (!fLowFsdClusters) {
        cout << "BmnDedx::Init(): branch BmnSiliconLowerCluster not found! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
 
    if (fRunId == -1) {
        cout << "BmnDedx::Init(): run ID not defined! Task will be deactivated" << endl;
        SetActive(kFALSE);
        return kERROR;
    }
 
    a = new Float_t**[fNSt];
    b = new Float_t**[fNSt];
 
    TString dir = getenv("VMCWORKDIR");
    TString calibFileName = dir + "/input/calibration.root";
    TFile* calibFile = TFile::Open(calibFileName);
 
    if (fVerbose > 1)
        cout << "Calibration file reading: " << calibFileName << endl;
 
    for (Int_t i = 0; i < fNSt; ++i) {
        a[i] = new Float_t*[fNSect[i]];
        b[i] = new Float_t*[fNSect[i]];
        for (Int_t j = 0; j < fNSect[i]; ++j) {
            a[i][j] = new Float_t[fNMomSlices];
            b[i][j] = new Float_t[fNMomSlices];
            for (Int_t k = 0; k < fNMomSlices; ++k) {
                TH1F* hA = (TH1F*)calibFile->Get(Form("hA_run_%d_%d_%d", i, j, k));
                TH1F* hB = (TH1F*)calibFile->Get(Form("hB_run_%d_%d_%d", i, j, k));
 
                a[i][j][k] = hA->GetBinContent(hA->FindBin(fRunId));
                b[i][j][k] = hB->GetBinContent(hB->FindBin(fRunId));
            }
        }
    }
    calibFile->Close();
 
    if (fVerbose > 1)
        cout << "\n======================== BmnDedx::Init finished ====================\n" << endl;
 
    return kSUCCESS;
}
 
void BmnDedx::CalcDedx(BmnGlobalTrack* tr)
{
    Float_t p = tr->GetP();
    if (p < 0) {
        tr->SetDedx(0.0);
        tr->SetQ(-1);
    }
    Int_t momIdx = -1;
    for (Int_t iP = 0; iP < fNMomSlices; ++iP) {
        if (p > fMinMomBound[iP] && p < fMaxMomBound[iP]) {
            momIdx = iP;
            break;
        }
    }
    if (momIdx == -1)
        return;
 
    CbmStsTrack* ststr = (CbmStsTrack*)fStsTracks->At(tr->GetGemTrackIndex());
    set<Float_t> signalsLow;
    Int_t nHits = ststr->GetNStsHits();
 
    Int_t nGemHits = 0;
    Int_t nFsdHits = 0;
    for (Int_t hitIdx = 0; hitIdx < nHits; ++hitIdx) {
        CbmStsHit* hit = (CbmStsHit*)fStsHits->At(ststr->GetStsHitIndex(hitIdx));
        Int_t sect = hit->GetSectorNr() - 1;
        Int_t st = hit->GetStationNr() - 1;
        StripCluster* low;
        if (hit->GetSystemId() == kGEM) {
            low = (StripCluster*)fLowGemClusters->At(hit->fDigiF % 1000000);
            nGemHits++;
        } else if (hit->GetSystemId() == kSILICON) {
            low = (StripCluster*)fLowFsdClusters->At(hit->fDigiF % 1000000);
            nFsdHits++;
        }
 
        signalsLow.insert(low->TotalSignal * a[st][sect][momIdx] + b[st][sect][momIdx]);
    }
 
    Int_t usedHits = (nHits == 3)    ? 2
                     : (nHits == 4)  ? 2
                     : (nHits == 5)  ? 3
                     : (nHits == 6)  ? 4
                     : (nHits == 7)  ? 4
                     : (nHits == 8)  ? 5
                     : (nHits == 9)  ? 6
                     : (nHits == 10) ? 6
                     : (nHits > 10)  ? 7
                                     : 0;
    auto it = signalsLow.begin();
    Float_t totSigLow = 0.0;
    for (Int_t i = 0; i < usedHits; ++i) {
        totSigLow += (*it);
        it++;
    }
 
    fDedx = (usedHits > 0) ? totSigLow / usedHits : 0.0;
 
    tr->SetDedx(fDedx);
    tr->SetQ(GetQ(tr));
 
    return;
}
 
Float_t BmnDedx::GetDedx(BmnGlobalTrack* tr)
{
    fDedx = tr->GetDedx();
    return fDedx;
}
 
Float_t BmnDedx::CalcNormalizedDedx(BmnGlobalTrack* tr)
{
    return fDedx / fAlephFit_He4->Eval(tr->GetP());
}
 
Short_t BmnDedx::GetQ(BmnGlobalTrack* tr)
{
    if (CalcNormalizedDedx(tr) > fNormDedxThresh)
        return 2;
    else
        return 1;
}
 
BmnDedx::~BmnDedx() {}