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 "UniDetectorParameter.h"
#include "uni_db_structures.h"
 
#include <TBufferJSON.h>
#include <TStopwatch.h>
#include <exception>
#include <iostream>
#include <nlohmann/json.hpp>
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)
{
    fnRun = 0;
    fnPeriod = 0;
 
    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");
 
        fDigiRunHeader = static_cast<DigiRunHeader*>(ioman->GetInFile()->Get("DigiRunHeader"));
        fnRun = fDigiRunHeader->GetRunId();
        fnPeriod = fDigiRunHeader->GetPeriodId();
        // std::cout << "fnRun = " << fnRun << "; fnPeriod = " << fnPeriod << std::endl;
        if (fnPeriod == 9) {
            fGr4FitBC1Top = new TGraph();
            fGr4FitBC1Bottom = new TGraph();
            f_gaus = new TF1("f_gaus", "gaus", 0, 500);
            fBC1Digits_Top = (TClonesArray*)ioman->GetObject("TQDC_BC1T");
            fBC1Digits_Bottom = (TClonesArray*)ioman->GetObject("TQDC_BC1B");
            GetT0SlewingFromDb();
        }
 
    } 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 (fnPeriod == 9 && isExpData && fBC1Digits_Top && fBC1Digits_Top->GetEntriesFast() > 0 && fBC1Digits_Bottom
            && fBC1Digits_Bottom->GetEntriesFast() > 0)
        {
            Float_t amp_temp, time_temp;
            time_temp = CalculateT0(&amp_temp);
            fTrigInfo->SetT0_time(time_temp);
            fTrigInfo->SetT0_amp(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();
    auto& valBC1T = digBC1T->GetValue();
    auto& valBC1B = digBC1B->GetValue();
    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.size() == 0 || valBC1B.size() == 0)
        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 = (int)LimitWvLow; ist < (int)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;
        }
    }
 
    if (maxBC1T < LimitLowAmpBC1T || maxBC1B < LimitLowAmpBC1B)
        return 0;
 
    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);
    }
}
 
// ---- clear points in TGraph ----------------------------------
bool BmnFillTrigInfoDst::GetT0SlewingFromDb()
{
    // defauilt parameters for run 9
    // The shape of the curve and its parameters depend on the properties of the T0 detectors
    f_BC1 = new TF1("f_BC1", "[0] + [1]*TMath::Exp([2]*([3]+x))", 0, 20000);
    f_BC1->SetParameters(0.5327, -3.611, -0.0002426, 855.5);
 
    // These constants depends on DAQ settings and WR sinchronization
    LimitTdcLow = 4225.;
    LimitTdcHgh = 4249.;
 
    // These constants depends on DAQ settings
    LimitWvLow = 320;
    LimitWvHgh = 350;
 
    // These constants depends on waveform of the signal
    GateFitLow = 5.5;
    GateFitHgh = 2.5;
 
    // These constants depends on waveform of the signal
    LimitLowAmpBC1T = 1000.;
    LimitLowAmpBC1B = 1000.;
 
    UniDetectorParameter* pDetectorParameter =
        UniDetectorParameter::GetDetectorParameter("T0", "T0_SlewingCorr", fnPeriod, fnRun);
    if (pDetectorParameter == nullptr) {
        LOG(info) << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t No data for T0 correction was found in the "
                     "database. The default values will be used.";
        return false;
    }
 
    JsonValue* json_value = (JsonValue*)pDetectorParameter->GetValue();
    Bool_t flagOK = true;
    // the record can be written by another ROOT version or be malformed, then the JSON conversion
    // throws an exception, which must not terminate the reconstruction: default values are used
    try {
        if (json_value->value.contains("SlewingFunction")) {
            // std::cout << json_value->value["SlewingFunction"] << std::endl;
            std::string s1 = json_value->value["SlewingFunction"].dump();
            TObject* slewingObject = TBufferJSON::ConvertFromJSON(s1.c_str());
            TF1* slewingFunction = dynamic_cast<TF1*>(slewingObject);
            if (slewingFunction != nullptr) {
                delete f_BC1;   // replace the default slewing function by the one from the database
                f_BC1 = slewingFunction;
            } else {
                delete slewingObject;
                LOG(error)
                    << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t The \"SlewingFunction\" field for T0 correction "
                       "was not read as TF1 from the database record. The default function will be used.";
                flagOK = false;
            }
        } else
            LOG(info)
                << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t The \"SlewingFunction\" field for T0 correction "
                   "was not found in the database record. The "
                   "default values will be used.";
 
        if (json_value->value.contains("Parameters")) {
            // std::cout << json_value->value["Parameters"] << std::endl;
 
            if (json_value->value["Parameters"].contains("LimitTdcLow"))
                LimitTdcLow = json_value->value["Parameters"]["LimitTdcLow"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("LimitTdcHgh"))
                LimitTdcHgh = json_value->value["Parameters"]["LimitTdcHgh"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("LimitWvLow"))
                LimitWvLow = json_value->value["Parameters"]["LimitWvLow"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("LimitWvHgh"))
                LimitWvHgh = json_value->value["Parameters"]["LimitWvHgh"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("GateFitLow"))
                GateFitLow = json_value->value["Parameters"]["GateFitLow"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("GateFitHgh"))
                GateFitHgh = json_value->value["Parameters"]["GateFitHgh"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("LimitLowAmpBC1T"))
                LimitLowAmpBC1T = json_value->value["Parameters"]["LimitLowAmpBC1T"];
            else
                flagOK = false;
 
            if (json_value->value["Parameters"].contains("LimitLowAmpBC1B"))
                LimitLowAmpBC1B = json_value->value["Parameters"]["LimitLowAmpBC1B"];
            else
                flagOK = false;
 
        } else
            LOG(info)
                << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t The \"Parameters\" field for T0 correction was "
                   "not found in the database record. The default "
                   "values will be used.";
    } catch (const std::exception& exc) {
        LOG(error) << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t The database record for T0 correction was not "
                      "parsed: "
                   << exc.what() << ". The default values will be used.";
        delete json_value;
        delete pDetectorParameter;
        return false;
    }
 
    if (flagOK == false)
        LOG(info) << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t Some parameters are not found in the "
                     "\"Parameters\" field. The default values will be used.";
 
    // clean memory after work
    delete json_value;
    delete pDetectorParameter;
 
    LOG(info) << "BmnFillTrigInfoDst::GetT0SlewingFromDb() :\n\t The parameters for calculating T0 have been "
                 "loaded from the database";
 
    return true;
}