BmnFillTrigInfoDst.cxx

Go to the documentation of this file.

/********************************************************************************
 *    BmnFillTrigInfoDst.cxx                                                    *
 *******************************************************************************/
 
#include "BmnFillTrigInfoDst.h"
 
#include "BmnFDPoint.h"
#include "BmnTrigDigit.h"
#include "BmnTrigInfoDst.h"
#include "BmnTrigWaveDigit.h"
#include "FairRootManager.h"
 
#include <TStopwatch.h>
#include <iostream>
using namespace std;
 
static Double_t workTime = 0.0;
 
// ---- Default constructor -------------------------------------------
BmnFillTrigInfoDst::BmnFillTrigInfoDst()
    : BmnFillTrigInfoDst(false)
{
    LOG(debug) << "Defaul Constructor of BmnFillTrigInfoDst";
}
 
// ---- Constructor with flag whether experimental data are used --------
BmnFillTrigInfoDst::BmnFillTrigInfoDst(Bool_t isExp)
    : BmnTask("BmnFillTrigInfoDst")
    , isExpData(isExp)
    , doProcessTriggers(true)
    , fDigiEventHeader(nullptr)
    , fVCDigits(nullptr)
    , fBC0Digits(nullptr)
    , fBC1Digits(nullptr)
    , fBC2Digits(nullptr)
    , fBDDigits(nullptr)
    , fSiMDDigits(nullptr)
    , fFDDigits(nullptr)
    , fBC1Digits_Top(nullptr)
    , fBC1Digits_Bottom(nullptr)
    , fTrigInfo(nullptr)
{
    fGr4FitBC1Top = new TGraph();
    fGr4FitBC1Bottom = new TGraph();
    f_gaus = new TF1("f_gaus", "gaus", 0, 500);
 
    // It is bad way to do it here. Need transmit run number to the task and get parameters from DB
    // These constants are derived from fitting the time-of-flight between BC2 and BC1.
    f_BC1 = new TF1("f_BC1", "[0] + [1]*TMath::Exp([2]*([3]+x))", 0, 20000);
    f_BC1->SetParameters(1.07489, -3.01848, -0.000169763, 358.637);
    // These constants depends on DAQ settings
    LimitTdcLow = 2600.;
    LimitTdcHgh = 2670.;
    LimitWvLow = 320;
    LimitWvHgh = 350;
    // These constants depends on waveform of the signal
    GateFitLow = 5.5;
    GateFitHgh = 2.5;
 
    LOG(debug) << "Constructor with experimental data flag of BmnFillTrigInfoDst";
}
 
// ---- Destructor ----------------------------------------------------
BmnFillTrigInfoDst::~BmnFillTrigInfoDst()
{
    LOG(debug) << "Destructor of BmnFillTrigInfoDst";
    if (fTrigInfo)
        delete fTrigInfo;
}
 
// ---- Init ----------------------------------------------------------
InitStatus BmnFillTrigInfoDst::Init()
{
    LOG(debug) << "Initialization of BmnFillTrigInfoDst";
 
    // Get a handle from the IO manager
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        LOG(error) << "Init: FairRootManager is not instantiated! BmnFillTrigInfoDst task will be deactivated";
        SetActive(kFALSE);
        return kERROR;
    }
 
    // Get input branches
    if (isExpData) {
        fDigiEventHeader = (BmnEventHeader*)ioman->GetObject("BmnEventHeader.");
        fVCDigits = (TClonesArray*)ioman->GetObject("TQDC_VCS");
        fBC0Digits = (TClonesArray*)ioman->GetObject("TQDC_BC0");
        fBC1Digits = (TClonesArray*)ioman->GetObject("TQDC_BC1S");
        if (!fBC1Digits)
            fBC1Digits = (TClonesArray*)ioman->GetObject("TQDC_BC1");
        fBC2Digits = (TClonesArray*)ioman->GetObject("TQDC_BC2AS");
        if (!fBC2Digits)
            fBC2Digits = (TClonesArray*)ioman->GetObject("TQDC_BC2");
        fBDDigits = (TClonesArray*)ioman->GetObject("BD");
        fSiMDDigits = (TClonesArray*)ioman->GetObject("SiMD");
        fFDDigits = (TClonesArray*)ioman->GetObject("TQDC_FD");
        fBC1Digits_Top = (TClonesArray*)ioman->GetObject("TQDC_BC1T");
        fBC1Digits_Bottom = (TClonesArray*)ioman->GetObject("TQDC_BC1B");
    } else {
        fBDDigits = (TClonesArray*)ioman->GetObject("BdDigit");
        fFDDigits = (TClonesArray*)ioman->GetObject("FDPoint");
    }
 
    if (isExpData && !fVCDigits && !fBC0Digits && !fBC1Digits && !fBC2Digits && !fBDDigits && !fSiMDDigits
        && !fFDDigits)
    {
        LOG(error) << "Init: all trigger branches not found! BmnFillTrigInfoDst task will process only basic data";
        doProcessTriggers = false;
        // SetActive(kFALSE);
        // return kERROR;
    }
 
    // Create and register TrigInfo object for trigger output data
    fTrigInfo = new BmnTrigInfoDst(isExpData);
    ioman->Register("BmnTrigInfo.", "TRIGGER", fTrigInfo, kTRUE);
    if (!fTrigInfo) {
        LOG(error) << "Init: BmnTrigInfo not created! BmnFillTrigInfoDst task will be deactivated";
        SetActive(kFALSE);
        return kERROR;
    }
 
    return kSUCCESS;
}
 
// ---- Exec ----------------------------------------------------------
void BmnFillTrigInfoDst::Exec(Option_t* /*option*/)
{
    if (!IsActive())
        return;
 
    TStopwatch sw;
    sw.Start();
 
    if (fVerbose > 1)
        cout << "\n======================== BmnFillTrigInfoDst exec started =====================\n" << endl;
 
    fTrigInfo->Reset();
 
    if (fDigiEventHeader) {
        fTrigInfo->SetInputSignalsAR(fDigiEventHeader->GetInputSignalsAR());
        fTrigInfo->SetInputSignalsBR(fDigiEventHeader->GetInputSignalsBR());
        fTrigInfo->SetInputSignalsVector(fDigiEventHeader->GetInputSignalsVector());
    }
 
    if (doProcessTriggers) {
        if (fVCDigits && fVCDigits->GetEntriesFast() > 0) {
            BmnTrigWaveDigit* digiVC = (BmnTrigWaveDigit*)fVCDigits->At(0);
            fTrigInfo->SetVCAmp(digiVC->GetPeak());
        }
        if (fBC0Digits && fBC0Digits->GetEntriesFast() > 0) {
            BmnTrigWaveDigit* digiBC0 = (BmnTrigWaveDigit*)fBC0Digits->At(0);
            fTrigInfo->SetBC0Amp(digiBC0->GetPeak());
            fTrigInfo->GetBC0Digits()->AbsorbObjects(fBC0Digits);
        }
        if (fBC1Digits && fBC1Digits->GetEntriesFast() > 0) {
            BmnTrigWaveDigit* digiBC1 = (BmnTrigWaveDigit*)fBC1Digits->At(0);
            fTrigInfo->SetBC1Integral(digiBC1->GetIntegral());
            fTrigInfo->SetBC1Amp(digiBC1->GetPeak());
            fTrigInfo->GetBC1Digits()->AbsorbObjects(fBC1Digits);
        }
        if (fBC2Digits && fBC2Digits->GetEntriesFast() > 0) {
            BmnTrigWaveDigit* digiBC2 = (BmnTrigWaveDigit*)fBC2Digits->At(0);
            fTrigInfo->SetBC2Amp(digiBC2->GetPeak());
            fTrigInfo->GetBC2Digits()->AbsorbObjects(fBC2Digits);
        }
        if (fBDDigits) {
            fTrigInfo->SetBDMult(fBDDigits->GetEntriesFast());
            fTrigInfo->GetBDDigits()->AbsorbObjects(fBDDigits);
        }
        if (fSiMDDigits)
            fTrigInfo->SetSiMDMult(fSiMDDigits->GetEntriesFast());
        if (fFDDigits && fFDDigits->GetEntriesFast() > 0) {
            if (isExpData) {
                BmnTrigWaveDigit* digiFD = (BmnTrigWaveDigit*)fFDDigits->At(0);
                fTrigInfo->SetFDAmp(digiFD->GetPeak());
            }
            fTrigInfo->GetFDDigits()->AbsorbObjects(fFDDigits);
        }
        if (isExpData && fBC1Digits_Top && fBC1Digits_Top->GetEntriesFast() > 0 && fBC1Digits_Bottom
            && fBC1Digits_Bottom->GetEntriesFast() > 0)
        {
            Float_t* amp_temp = new Float_t();
            fTrigInfo->SetT0_time(CalculateT0(amp_temp));
            fTrigInfo->SetT0_amp(*amp_temp);
            delete amp_temp;
        } else {
            fTrigInfo->SetT0_time(0.);
            fTrigInfo->SetT0_amp(0.);
        }
    }
 
    sw.Stop();
    workTime += sw.RealTime();
}
 
// ---- Finish --------------------------------------------------------
void BmnFillTrigInfoDst::Finish()
{
    LOG(debug) << "Finish of BmnFillTrigInfoDst";
    printf("Work time of BmnFillTrigInfoDst: %4.2f sec.\n", workTime);
}
 
// ---- CalculateT0 ----------------------------------------------------
Float_t BmnFillTrigInfoDst::CalculateT0(Float_t* amp)
{
    Float_t tempT0 = 0.;
    if (amp)
        *amp = 0;
 
    // Get WaveDigit for BC1_Top and BC1_Bottom detectors
    BmnTrigWaveDigit* digBC1T = (BmnTrigWaveDigit*)fBC1Digits_Top->UncheckedAt(0);
    BmnTrigWaveDigit* digBC1B = (BmnTrigWaveDigit*)fBC1Digits_Bottom->UncheckedAt(0);
    Int_t nSemplesWF = digBC1T->GetNSamples();
    Short_t* valBC1T = digBC1T->GetShortValue();
    Short_t* valBC1B = digBC1B->GetShortValue();
    vector<Double_t> vtimeBC1T = digBC1T->TdcVector();
    vector<Double_t> vtimeBC1B = digBC1B->TdcVector();
 
    if (nSemplesWF <= 0 || LimitWvHgh >= nSemplesWF)
        return 0;
 
    if (vtimeBC1T.size() == 0 || vtimeBC1B.size() == 0)
        return 0;
 
    if (valBC1T == nullptr || valBC1B == nullptr)
        return 0;
 
    Double_t timeBC1 = 0, timeBC1T = 0, timeBC1B = 0;
 
    // Find TDC time for Top and Bottom in TdcDigitsArray
    timeBC1T = FindTimeFromVector(vtimeBC1T);
    timeBC1B = FindTimeFromVector(vtimeBC1B);
    timeBC1 = (timeBC1T + timeBC1B) * 0.5;
 
    // calculate the signal amplitude of Top and Bottom
    Double_t maxBC1T = 0, maxBC1B = 0, maxBC1T_fit = 0, maxBC1B_fit = 0;
    Double_t maxStBC1T = 0, maxStBC1B = 0;
 
    ClearTGraph(fGr4FitBC1Top);
    ClearTGraph(fGr4FitBC1Bottom);
    for (Int_t ist = LimitWvLow; ist < LimitWvHgh; ist++) {
        fGr4FitBC1Top->AddPoint(ist, valBC1T[ist]);
        fGr4FitBC1Bottom->AddPoint(ist, valBC1B[ist]);
        if (maxBC1T < valBC1T[ist]) {
            maxBC1T = valBC1T[ist];
            maxStBC1T = (Double_t)ist;
        }
        if (maxBC1B < valBC1B[ist]) {
            maxBC1B = valBC1B[ist];
            maxStBC1B = (Double_t)ist;
        }
    }
 
    f_gaus->SetParameters(maxBC1T, maxStBC1T, 1);
    fGr4FitBC1Top->Fit(f_gaus, "QW", "", maxStBC1T - GateFitLow, maxStBC1T + GateFitHgh);
    maxBC1T_fit = f_gaus->GetParameter(0);
 
    f_gaus->SetParameters(maxBC1B, maxStBC1B, 1);
    fGr4FitBC1Bottom->Fit(f_gaus, "QW", "", maxStBC1B - GateFitLow, maxStBC1B + GateFitHgh);
    maxBC1B_fit = f_gaus->GetParameter(0);
 
    // calculate the geometric mean of amplitude
    Double_t gAmpBC1 = TMath::Sqrt(maxBC1T_fit * maxBC1B_fit);
 
    // correct the time on amplitude
    // tempT0 = timeBC1;
    tempT0 = timeBC1 + f_BC1->Eval(gAmpBC1);
    if (amp)
        *amp = gAmpBC1;
 
    return tempT0;
}
 
// ---- find time from vector of TDC hits ----------------------------------
Double_t BmnFillTrigInfoDst::FindTimeFromVector(vector<Double_t> vtime)
{
 
    for (UInt_t i = 0; i < vtime.size(); i++) {
        if (vtime.at(i) > LimitTdcLow && vtime.at(i) < LimitTdcHgh)
            return vtime.at(i);
    }
    return 0;
}
 
// ---- clear points in TGraph ----------------------------------
void BmnFillTrigInfoDst::ClearTGraph(TGraph* g1)
{
    if (g1 == nullptr)
        return;
 
    for (int ip = g1->GetN() - 1; ip >= 0; ip--) {
        g1->RemovePoint(ip);
    }
}