BmnSiResiduals.cxx

Go to the documentation of this file.

#include "BmnSiResiduals.h"
 
BmnSiResiduals::BmnSiResiduals(Int_t run_period, Int_t run_number, Double_t fieldScale)
    : isField(kFALSE)
    , fDebug(kFALSE)
    , outRes(NULL)
    , isPrintToFile(kFALSE)
    , isMergedDigits(kFALSE)
{
    fPeriod = run_period;
    fNumber = run_number;
 
    if (Abs(fieldScale) > DBL_EPSILON)
        isField = kTRUE;
 
    TString gPathSiliconConfig = gSystem->Getenv("VMCWORKDIR");
    gPathSiliconConfig += "/parameters/silicon/XMLConfigs/";
 
    // Create Silicon detector --------------------------------------------------
    TString config = (fPeriod == 6) ? "SiliconRunSpring2017.xml" : (fPeriod == 7) ? "SiliconRunSpring2018.xml" : "";
 
    fDetector = new BmnSiliconStationSet(gPathSiliconConfig + config);
 
    fBranchSiHits = "BmnSiliconHit";
    fBranchGlobalTracks = "BmnGlobalTrack";
    fBranchResiduals = "BmnResiduals";
    fBranchFairEventHeader = "EventHeader.";
}
 
void BmnSiResiduals::Finish()
{
    delete fDetector;
}
 
InitStatus BmnSiResiduals::Init()
{
    FairRootManager* ioman = FairRootManager::Instance();
 
    fSiHits = (TClonesArray*)ioman->GetObject(fBranchSiHits.Data());
    fGlobalTracks = (TClonesArray*)ioman->GetObject(fBranchGlobalTracks.Data());
    fFairEventHeader = (FairEventHeader*)ioman->GetObject(fBranchFairEventHeader.Data());
 
    fSiResiduals = new TClonesArray(fBranchResiduals.Data());
 
    ioman->Register("BmnSiResiduals", "RESID", fSiResiduals, kTRUE);
    if (TString(ioman->GetInFile()->GetName()).Contains("merge"))
        isMergedDigits = kTRUE;
}
 
void BmnSiResiduals::Exec(Option_t* opt)
{
    //  fFairEventHeader->SetRunId((isMergedDigits) ? 0 : fNumber);
    fSiResiduals->Delete();
    Residuals();
}
 
void BmnSiResiduals::Residuals()
{
    for (Int_t iGlobTrack = 0; iGlobTrack < fGlobalTracks->GetEntriesFast(); iGlobTrack++) {
        BmnGlobalTrack* glTrack = (BmnGlobalTrack*)fGlobalTracks->UncheckedAt(iGlobTrack);
 
        Double_t Tx = glTrack->GetParamFirst()->GetTx();
        Double_t X0 = glTrack->GetParamFirst()->GetX();
        Double_t Ty = glTrack->GetParamFirst()->GetTy();
        Double_t Y0 = glTrack->GetParamFirst()->GetY();
 
        Double_t Z0 = glTrack->GetParamFirst()->GetZ();
 
        //        Int_t silIdx = glTrack->GetSilHitIndex();
        //        if (silIdx == -1)
        //            return;
        vector<Int_t> silIndices = glTrack->GetSilHitIndices();
        if (silIndices.size() == 0)
            return;
 
        for (Int_t iHitIdx = 0; iHitIdx < silIndices.size(); iHitIdx++) {
            Double_t xRes = 0., yRes = 0.;
            BmnSiliconHit* hit = (BmnSiliconHit*)fSiHits->UncheckedAt(silIndices.at(iHitIdx));
            Double_t x = hit->GetX();
            Double_t y = hit->GetY();
            Double_t z = hit->GetZ();
 
            Int_t stat = hit->GetStation();
            Int_t mod = hit->GetModule();
 
            if (!isField) {
                xRes = x - (X0 + Tx * (z - Z0));
                yRes = y - (Y0 + Ty * (z - Z0));
            }
 
            BmnResiduals* resid =
                new ((*fSiResiduals)[fSiResiduals->GetEntriesFast()]) BmnResiduals(stat, mod, xRes, yRes, 0., isField);
            resid->SetTrackId(iGlobTrack);
            resid->SetHitId(silIndices.at(iHitIdx));
        }
    }
}