TRestAnalysisTree.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
 
#include "TRestAnalysisTree.h"
 
#include <TFile.h>
#include <TH1F.h>
#include <TLeaf.h>
#include <TObjArray.h>
 
#include "TRestStringHelper.h"
#include "TRestStringOutput.h"
 
using namespace std;
 
ClassImp(TRestAnalysisTree) using namespace std;
 
TRestAnalysisTree::TRestAnalysisTree() : TTree() {
    SetName("AnalysisTree");
    SetTitle("REST Analysis tree");
 
    Initialize();
}
 
TRestAnalysisTree::TRestAnalysisTree(TTree* tree) : TTree() {
    SetName("AnalysisTree");
    SetTitle("Linked to " + (TString)tree->GetName());
 
    Initialize();
 
    tree->GetEntry(0);
 
    TObjArray* lfs = tree->GetListOfLeaves();
    for (int i = 0; i < lfs->GetLast() + 1; i++) {
        TLeaf* lf = (TLeaf*)lfs->UncheckedAt(i);
        RESTValue obs;
        ReadLeafValueToObservable(lf, obs);
        obs.name = lf->GetName();
        SetObservable(-1, obs);
    }
}
 
TRestAnalysisTree* TRestAnalysisTree::ConvertFromTTree(TTree* tree) {
    if (tree->ClassName() == (TString) "TRestAnalysisTree") {
        return (TRestAnalysisTree*)tree;
    }
 
    return new TRestAnalysisTree(tree);
}
 
TRestAnalysisTree::TRestAnalysisTree(const TString& name, const TString& title) : TTree(name, title) {
    SetName(name);
    SetTitle(title);
 
    Initialize();
}
 
void TRestAnalysisTree::Initialize() {
    fRunOrigin = 0;
    fSubRunOrigin = 0;
    fEventID = 0;
    fSubEventID = 0;
    fSubEventTag = new TString();
    fTimeStamp = 0;
 
    fStatus = 0;
    fNObservables = 0;
    fChain = nullptr;
 
    fObservableIdMap.clear();
    fObservableIdSearchMap.clear();
    fObservableDescriptions.clear();
    fObservableNames.clear();
    fObservables.clear();
}
 
Int_t TRestAnalysisTree::GetObservableID(const string& obsName) {
    MakeObservableIdMap();
    auto iter = fObservableIdMap.find(obsName);
    if (iter != fObservableIdMap.end()) {
        return iter->second;
    } else {
        return -1;
    }
}
 
Int_t TRestAnalysisTree::GetMatchedObservableID(const string& obsName) {
    // if (ObservableExists(obsName)) return GetObservableID(obsName);
    auto iter = fObservableIdSearchMap.find(obsName);
    if (iter != fObservableIdSearchMap.end()) {
        return iter->second;
    } else {
        vector<int> diffs;
        for (auto& fObservableName : fObservableNames) {
            string obs = (string)fObservableName;
            int diff1 = DiffString(obs, obsName);
 
            obs = obs.substr(obs.find('_') + 1, -1);
            int diff2 = DiffString(obs, obsName);
 
            diffs.push_back(min(diff1, diff2));
        }
 
        int matchedCount = 0;
        int minPos = 0;
        int min = 1e5;
        for (unsigned int i = 0; i < diffs.size(); i++) {
            if (min > diffs[i]) {
                min = diffs[i];
                minPos = i;
            }
            if (diffs[i] == 0) {
                matchedCount++;
            }
        }
 
        if (matchedCount == 1) {
            cout << "TRestAnalysisTree::GetObservableID(): Find observable " << fObservableNames[minPos]
                 << " for obs name: " << obsName << endl;
            fObservableIdSearchMap[obsName] = minPos;
            return minPos;
        } else if (matchedCount == 0) {
            RESTError << "TRestAnalysisTree::GetObservableID(): No Matching observable found" << RESTendl;
            if (min <= 2) {
                RESTError << "did you mean \"" << fObservableNames[minPos] << "\" ?" << RESTendl;
                fObservableIdSearchMap[obsName] = -1;
                return -1;
            }
        } else {
            RESTWarning << "TRestAnalysisTree::GetObservableID(): Multiple observables found" << RESTendl;
            RESTWarning
                << "Please specify the full observable name (e.g, sAna_ThresholdIntegral) in your code. "
                << RESTendl;
            RESTWarning << "Returning the first matched observable : " << fObservableNames[minPos]
                        << RESTendl;
            fObservableIdSearchMap[obsName] = minPos;
            return minPos;
        }
    }
    fObservableIdSearchMap[obsName] = -1;
    return -1;
}
 
Bool_t TRestAnalysisTree::ObservableExists(const string& obsName) {
    MakeObservableIdMap();
    return fObservableIdMap.count(obsName) > 0;
}
 
// Status              | Nbranches  |  NObservables  | Entries |  Description
// 1. Created          |     0      |       0        |    0    |  The AnalysisTree is just created in code
// 2. Retrieved        |    >=6     |       0        |    >0   |  The AnalysisTree is just retrieved from file
// 3. EmptyCloned      |    >=6     |       0        |    0    |  The AnalysisTree is cloned from another
// tree, it is empty
// 4. Connected        |    >=6     |  >=Nbranches-6 |    0    |  Local observables are created and connected
// to the branches. We can still modify observables in this case.
// 5. Filled           |    >=6     |     =NObs      |    >0   |  Local observables are created and connected
// to the branches, and entry > 0. We cannot modify observables in this case. if Nbranches == 6, We regard the
// tree as Connected/Filled
//
// Status              |  AddObs Logic  |  Fill Logic  | GetEntry Logic
// 1. Created          |    Allowed     |   ->4, ->5   |    ->4
// 2. Retrieved        |    Denied      |     ->5      |    ->5
// 3. EmptyCloned      |     ->4        |   ->4, ->5   |    ->4
// 4. Connected        |    Allowed     |     ->5      |    Allowed
// 5. Filled           |    Denied      |     Allowed  |    Allowed
//
int TRestAnalysisTree::EvaluateStatus() {
    int NBranches = GetListOfBranches()->GetEntriesFast();
    int NObsInList = fObservables.size();
    int Entries = GetEntriesFast();
 
    // if (fROOTTree != nullptr) {
    //    return TRestAnalysisTree_Status::ROOTTree;
    //}
 
    if (NBranches == 0) {
        return TRestAnalysisTree_Status::Created;
    }
 
    if (NBranches > 6) {
        if (NObsInList == 0) {
            if (Entries > 0) {
                return TRestAnalysisTree_Status::Retrieved;
            } else {
                return TRestAnalysisTree_Status::EmptyCloned;
            }
        } else {
            if (Entries > 0) {
                return TRestAnalysisTree_Status::Filled;
            } else {
                return TRestAnalysisTree_Status::Connected;
            }
        }
    } else if (NBranches == 6) {
        auto br = (TBranch*)GetListOfBranches()->UncheckedAt(0);
        if ((int*)br->GetAddress() != &fRunOrigin) {
            if (Entries > 0) {
                return TRestAnalysisTree_Status::Retrieved;
            } else {
                return TRestAnalysisTree_Status::EmptyCloned;
            }
        } else {
            if (Entries > 0) {
                return TRestAnalysisTree_Status::Filled;
            } else {
                return TRestAnalysisTree_Status::Connected;
            }
        }
    } else {
        return TRestAnalysisTree_Status::Error;
    }
 
    return TRestAnalysisTree_Status::Error;
}
 
void TRestAnalysisTree::UpdateObservables() {
    // connect basic event branches
    TBranch* br1 = GetBranch("runOrigin");
    TBranch* br2 = GetBranch("subRunOrigin");
    TBranch* br3 = GetBranch("eventID");
    TBranch* br4 = GetBranch("subEventID");
    TBranch* br5 = GetBranch("timeStamp");
    TBranch* br6 = GetBranch("subEventTag");
 
    if (br1 && br2 && br3 && br4 && br5 && br6) {
        br1->SetAddress(&fRunOrigin);
        br2->SetAddress(&fSubRunOrigin);
        br3->SetAddress(&fEventID);
        br4->SetAddress(&fSubEventID);
        br5->SetAddress(&fTimeStamp);
        br6->SetAddress(&fSubEventTag);
    } else {
        cout << "REST Error: TRestAnalysisTree::ConnectBranches(): event branches does not exist!" << endl;
        exit(1);
    }
 
    // create observables(no assembly) from stored observable info.
    fObservables = std::vector<RESTValue>(GetNumberOfObservables());
    for (int i = 0; i < GetNumberOfObservables(); i++) {
        fObservables[i] = REST_Reflection::WrapType((string)fObservableTypes[i]);
        fObservables[i].name = fObservableNames[i];
    }
    MakeObservableIdMap();
 
    if (fChain) {
        fChain->GetEntry(0);
    } else {
        TTree::GetEntry(0);
    }
 
    for (int i = 0; i < GetNumberOfObservables(); i++) {
        TBranch* branch = GetBranch(fObservableNames[i]);
        if (branch != nullptr) {
            if (branch->GetAddress() != nullptr) {
                if ((string)branch->ClassName() != "TBranch") {
                    // for TBranchElement the saved address is char**
                    fObservables[i].address = *(char**)branch->GetAddress();
                } else {
                    // for TBranch the saved address is char*
                    fObservables[i].address = branch->GetAddress();
                }
            } else {
                fObservables[i].Assembly();
                branch->SetAddress(fObservables[i].address);
            }
        }
    }
 
    fStatus = EvaluateStatus();
}
 
void TRestAnalysisTree::UpdateBranches() {
    if (!GetBranch("runOrigin")) Branch("runOrigin", &fRunOrigin);
    if (!GetBranch("subRunOrigin")) Branch("subRunOrigin", &fSubRunOrigin);
    if (!GetBranch("eventID")) Branch("eventID", &fEventID);
    if (!GetBranch("subEventID")) Branch("subEventID", &fSubEventID);
    if (!GetBranch("timeStamp")) Branch("timeStamp", &fTimeStamp);
    if (!GetBranch("subEventTag")) Branch("subEventTag", fSubEventTag);
 
    for (int n = 0; n < GetNumberOfObservables(); n++) {
        TString typeName = fObservableTypes[n];
        TString brName = fObservableNames[n];
        char* ref = fObservables[n].address;
 
        TBranch* branch = GetBranch(brName);
        if (branch == nullptr) {
            if (typeName == "double") {
                this->Branch(brName, (double*)ref);
            } else if (typeName == "float") {
                this->Branch(brName, (float*)ref);
            } else if (typeName == "long double") {
                this->Branch(brName, (long double*)ref);
            } else if (typeName == "bool") {
                this->Branch(brName, (bool*)ref);
            } else if (typeName == "char") {
                this->Branch(brName, (char*)ref);
            } else if (typeName == "int") {
                this->Branch(brName, (int*)ref);
            } else if (typeName == "short") {
                this->Branch(brName, (short*)ref);
            } else if (typeName == "unsigned char") {
                this->Branch(brName, (unsigned char*)ref);
            } else if (typeName == "unsigned int") {
                this->Branch(brName, (unsigned int*)ref);
            } else if (typeName == "unsigned short") {
                this->Branch(brName, (unsigned short*)ref);
            } else if (typeName == "long") {
                this->Branch(brName, (long*)ref);
            } else if (typeName == "long long") {
                this->Branch(brName, (long long*)ref);
            } else if (typeName == "unsigned long") {
                this->Branch(brName, (unsigned long*)ref);
            } else if (typeName == "unsigned long long") {
                this->Branch(brName, (unsigned long long*)ref);
            } else {
                this->Branch(brName, typeName, ref);
            }
 
            this->GetBranch(brName)->SetTitle("(" + typeName + ") " + fObservableDescriptions[n]);
        }
    }
 
    fStatus = EvaluateStatus();
}
 
void TRestAnalysisTree::InitObservables() {
    fObservables = std::vector<RESTValue>(GetNumberOfObservables());
    for (int i = 0; i < GetNumberOfObservables(); i++) {
        fObservables[i] = REST_Reflection::WrapType((string)fObservableTypes[i]);
        fObservables[i].name = fObservableNames[i];
    }
    MakeObservableIdMap();
}
 
void TRestAnalysisTree::MakeObservableIdMap() {
    if (fObservableIdMap.size() != fObservableNames.size()) {
        fObservableIdMap.clear();
        for (unsigned int i = 0; i < fObservableNames.size(); i++) {
            fObservableIdMap[(string)fObservableNames[i]] = i;
        }
        if (fObservableIdMap.size() != fObservableNames.size()) {
            cout << "REST ERROR! duplicated or blank observable name!" << endl;
        }
    }
}
 
void TRestAnalysisTree::ReadLeafValueToObservable(TLeaf* lf, RESTValue& obs) {
    if (lf == nullptr || lf->GetLen() == 0) return;
 
    if (obs.IsZombie()) {
        // this means we need to determine the observable type first
        string type;
        string leafdef(lf->GetBranch()->GetTitle());
        if (leafdef.find("[") == string::npos) {
            // there is no [] mark in the leaf. The leaf is a single values
            switch (leafdef[leafdef.size() - 1]) {
                case 'C':
                    type = "string";
                    break;
                case 'B':
                    type = "char";
                    break;
                case 'S':
                    type = "stort";
                    break;
                case 'I':
                    type = "int";
                    break;
                case 'F':
                    type = "float";
                    break;
                case 'D':
                    type = "double";
                    break;
                case 'L':
                    type = "long long";
                    break;
                case 'O':
                    type = "bool";
                    break;
                default:
                    RESTError << "Unsupported leaf definition: " << leafdef << "! Assuming int" << RESTendl;
                    type = "int";
            }
        } else {
            switch (leafdef[leafdef.size() - 1]) {
                case 'I':
                    type = "vector<int>";
                    break;
                case 'F':
                    type = "vector<float>";
                    break;
                case 'D':
                    type = "vector<double>";
                    break;
                default:
                    RESTError << "Unsupported leaf definition: " << leafdef << "! Assuming vector<int>"
                              << RESTendl;
                    type = "vector<int>";
            }
        }
        obs = REST_Reflection::Assembly(type);
    }
 
    // start to fill value
    if (obs.is_data_type) {
        // the observable is basic type, we directly set it address from leaf
        obs.address = (char*)lf->GetValuePointer();
    } else if (obs.type == "vector<double>") {
        auto ptr = (double*)lf->GetValuePointer();
        vector<double> val(ptr, ptr + lf->GetLen());
        obs.SetValue(val);
    } else if (obs.type == "vector<int>") {
        int* ptr = (int*)lf->GetValuePointer();
        vector<int> val(ptr, ptr + lf->GetLen());
        obs.SetValue(val);
    } else if (obs.type == "vector<float>") {
        auto ptr = (float*)lf->GetValuePointer();
        vector<float> val(ptr, ptr + lf->GetLen());
        obs.SetValue(val);
    } else if (obs.type == "string") {
        char* ptr = (char*)lf->GetValuePointer();
        string val(ptr, lf->GetLen());
        obs.SetValue(val);
    } else {
        RESTWarning << "Unsupported observable type to convert from TLeaf!" << RESTendl;
    }
}
 
RESTValue TRestAnalysisTree::GetObservable(const string& obsName) {
    Int_t id = GetObservableID(obsName);
    if (id == -1) return {};
    return GetObservable(id);
}
 
RESTValue TRestAnalysisTree::GetObservable(Int_t n) {
    if (n >= fNObservables) {
        cout << "Error! TRestAnalysisTree::GetObservable(): index outside limits!" << endl;
        return {};
    }
 
    if (fChain != nullptr) {
        return ((TRestAnalysisTree*)fChain->GetTree())->GetObservable(n);
    }
 
    return fObservables[n];
}
TString TRestAnalysisTree::GetObservableName(Int_t n) {
    if (n >= fNObservables) {
        cout << "Error! TRestAnalysisTree::GetObservableName(): index outside limits!" << endl;
        return "";
    }
    return fObservableNames[n];
}
TString TRestAnalysisTree::GetObservableDescription(Int_t n) {
    if (n >= fNObservables) {
        cout << "Error! TRestAnalysisTree::GetObservableDescription(): index outside limits!" << endl;
        return "";
    }
    return fObservableDescriptions[n];
}
TString TRestAnalysisTree::GetObservableType(Int_t n) {
    if (n >= fNObservables) {
        cout << "Error! TRestAnalysisTree::GetObservableType(): index outside limits!" << endl;
        return "double";
    }
    return fObservableTypes[n];
}
TString TRestAnalysisTree::GetObservableType(const string& obsName) {
    Int_t id = GetObservableID(obsName);
    if (id == -1) return "NotFound";
    return GetObservableType(id);
}
Double_t TRestAnalysisTree::GetDblObservableValue(const string& obsName) {
    return GetDblObservableValue(GetObservableID(obsName));
}
Double_t TRestAnalysisTree::GetDblObservableValue(Int_t n) {
    if (n >= fNObservables) {
        cout << "Error! TRestAnalysisTree::GetDblObservableValue(): index outside limits!" << endl;
        return 0;
    }
 
    if (GetObservableType(n) == "int") return GetObservableValue<int>(n);
    if (GetObservableType(n) == "float") return GetObservableValue<float>(n);
    if (GetObservableType(n) == "double") return GetObservableValue<double>(n);
 
    cout << "TRestAnalysisTree::GetDblObservableValue. Type " << GetObservableType(n)
         << " not supported! Returning zero" << endl;
    return 0.;
}
 
RESTValue TRestAnalysisTree::AddObservable(const TString& observableName, const TString& observableType,
                                           const TString& description) {
    if (fStatus == None) fStatus = EvaluateStatus();
 
    if (fStatus == Created) {
    } else if (fStatus == Retrieved) {
        cout << "REST_Warning: cannot add observable for retrieved tree with data!" << endl;
        return -1;
    } else if (fStatus == EmptyCloned) {
        UpdateObservables();
    } else if (fStatus == Connected) {
    } else if (fStatus == Filled) {
        cout << "REST_Warning: cannot add observable! AnalysisTree is already filled" << endl;
        return -1;
    } else if (fChain != nullptr) {
        cout << "Error! cannot add observable! AnalysisTree is in chain state" << endl;
        return -1;
    }
 
    if (GetObservableID((string)observableName) == -1) {
        RESTValue ptr = REST_Reflection::Assembly((string)observableType);
        ptr.name = observableName;
        if (!ptr.IsZombie()) {
            fObservableNames.push_back(observableName);
            fObservableIdMap[(string)observableName] = fObservableNames.size() - 1;
            fObservableDescriptions.push_back(description);
            fObservableTypes.push_back(observableType);
            fObservables.push_back(ptr);
 
            fNObservables++;
        } else {
            cout << "Error adding observable \"" << observableName << "\" with type \"" << observableType
                 << "\", type not found!" << endl;
            return -1;
        }
    } else {
        return GetObservableID((string)observableName);
    }
 
    return fObservables[fNObservables - 1];
}
 
void TRestAnalysisTree::PrintObservables() {
    if (fStatus == None) fStatus = EvaluateStatus();
 
    if (fStatus == Created) {
    } else if (fStatus == Retrieved) {
        cout << "TRestAnalysisTree::PrintObservables(): Reading entry 0" << endl;
        GetEntry(0);
    } else if (fStatus == EmptyCloned) {
        cout << "REST_Warning: AnalysisTree is empty" << endl;
        UpdateObservables();
    } else if (fStatus == Connected) {
        cout << "REST_Warning: AnalysisTree is empty" << endl;
    } else if (fStatus == Filled) {
    }
 
    cout.precision(15);
    std::cout << "Entry : " << GetReadEntry() << std::endl;
    std::cout << "> Event ID : " << GetEventID() << std::endl;
    std::cout << "> Event Time : " << GetTimeStamp() << std::endl;
    std::cout << "> Event Tag : " << GetSubEventTag() << std::endl;
    std::cout << "-----------------------------------------" << std::endl;
    cout.precision(6);
 
    for (int n = 0; n < GetNumberOfObservables(); n++) {
        PrintObservable(n);
    }
}
 
// print the Nth observable in the observal list
void TRestAnalysisTree::PrintObservable(int n) {
    if (n < 0 || n >= fNObservables) {
        return;
    }
    string obsVal = GetObservable(n).ToString();
    int lengthRemaining = Console::GetWidth() - 14 - 45 - 13;
 
    std::cout << "Observable : " << ToString(fObservableNames[n], 45)
              << "    Value : " << ToString(obsVal, lengthRemaining) << std::endl;
}
 
Int_t TRestAnalysisTree::GetEntry(Long64_t entry, Int_t getall) {
    if (fStatus == None) fStatus = EvaluateStatus();
 
    if (fStatus == Created) {
        UpdateBranches();
    } else if (fStatus == Retrieved) {
        UpdateObservables();
    } else if (fStatus == EmptyCloned) {
        UpdateObservables();
    } else if (fStatus == Connected) {
    } else if (fStatus == Filled) {
    }
 
    if (fChain != nullptr) {
        return fChain->GetEntry(entry, getall);
    } else {
        return TTree::GetEntry(entry, getall);
    }
}
 
void TRestAnalysisTree::SetEventInfo(TRestAnalysisTree* tree) {
    if (fChain != nullptr) {
        cout << "Error! cannot fill tree. AnalysisTree is in chain state" << endl;
        return;
    }
 
    if (tree != nullptr) {
        fEventID = tree->GetEventID();
        fSubEventID = tree->GetSubEventID();
        fTimeStamp = tree->GetTimeStamp();
        *fSubEventTag = tree->GetSubEventTag();
        fRunOrigin = tree->GetRunOrigin();
        fSubRunOrigin = tree->GetSubRunOrigin();
    }
}
 
void TRestAnalysisTree::SetEventInfo(TRestEvent* evt) {
    if (fChain != nullptr) {
        cout << "Error! cannot fill tree. AnalysisTree is in chain state" << endl;
        return;
    }
 
    if (evt != nullptr) {
        fEventID = evt->GetID();
        fSubEventID = evt->GetSubID();
        fTimeStamp = evt->GetTimeStamp().AsDouble();
        *fSubEventTag = evt->GetSubEventTag();
        fRunOrigin = evt->GetRunOrigin();
        fSubRunOrigin = evt->GetSubRunOrigin();
    }
}
 
Int_t TRestAnalysisTree::Fill() {
    if (fStatus == None) fStatus = EvaluateStatus();
 
    if (fChain != nullptr) {
        cout << "Error! cannot fill tree. AnalysisTree is in chain state" << endl;
        return -1;
    }
 
    fSetObservableIndex = 0;
 
    if (fStatus == Filled) {
        return TTree::Fill();
    } else if (fStatus == Created) {
        UpdateBranches();
        int a = TTree::Fill();
        fStatus = EvaluateStatus();
        return a;
    } else if (fStatus == Retrieved) {
        UpdateObservables();
        int a = TTree::Fill();
        fStatus = EvaluateStatus();
        return a;
    } else if (fStatus == EmptyCloned) {
        UpdateObservables();
        int a = TTree::Fill();
        fStatus = EvaluateStatus();
        return a;
    } else if (fStatus == Connected) {
        UpdateBranches();
        int a = TTree::Fill();
        fStatus = EvaluateStatus();
        return a;
    }
    return -1;
}
 
void TRestAnalysisTree::SetObservable(Int_t id, RESTValue obs) {
    if (fChain != nullptr) {
        cout << "Error! cannot set observable! AnalysisTree is in chain state" << endl;
        return;
    }
    if (id == -1) {
        // this means we want to find observable id by its name
        id = GetObservableID(obs.name);
        // if not found, we have a chance to create observable
        if (id == -1) {
            AddObservable(obs.name, obs.type);
            id = GetObservableID(obs.name);
        }
    } else if (id == (int)fObservables.size()) {
        // this means we want to add observable directly
        AddObservable(obs.name, obs.type);
        id = GetObservableID(obs.name);
    }
    if (id != -1) {
        if (obs.type != fObservables[id].type) {
            // if the observable branches are not created, and the type doesn't match,
            // we still have the chance to fix. We reset fObservableTypes and fObservableMemory
            // according to the input type value.
            cout << "Warning: SetObservableValue(): setting value with different type. Observable : "
                 << obs.name << ", existing type: " << fObservables[id].type
                 << ", value to add is in type: " << obs.type << ", trying to update to the new type."
                 << endl;
            fObservableTypes[id] = obs.type;
            string name = fObservables[id].name;
            fObservables[id].Destroy();
            fObservables[id] = REST_Reflection::Assembly(obs.type);
            fObservables[id].name = name;
        }
        obs >> fObservables[id];
    }
}
 
void TRestAnalysisTree::SetObservable(string name, RESTValue value) {
    value.name = name;
    SetObservable(-1, value);
}
 
Bool_t TRestAnalysisTree::EvaluateCuts(const string& cut) {
    std::vector<string> cuts = Split(cut, "&&", false, true);
 
    for (auto& cut : cuts) {
        if (!EvaluateCut(cut)) {
            return false;
        }
    }
 
    return true;
}
 
Bool_t TRestAnalysisTree::EvaluateCut(const string& cut) {
    const std::vector<string> validOperators = {"==", "<=", ">=", "!=", "=", ">", "<"};
 
    if (cut.find("&&") != string::npos) return EvaluateCuts(cut);
 
    string operation, observable;
    Double_t value = -1;
    for (const auto& validOperator : validOperators) {
        if (cut.find(validOperator) != string::npos) {
            operation = validOperator;
            observable = (string)cut.substr(0, cut.find(operation));
            value = std::stod((string)cut.substr(cut.find(operation) + operation.length(), string::npos));
            break;
        }
    }
 
    if (operation.empty())
        cout << "TRestAnalysisTree::EvaluateCut. Invalid operator in cut definition! " << cut << endl;
 
    Double_t val = GetDblObservableValue(observable);
    if (operation == "==" && value == val) return true;
    if (operation == "!=" && value != val) return true;
    if (operation == "<=" && val <= value) return true;
    if (operation == ">=" && val >= value) return true;
    if (operation == "=" && val == value) return true;
    if (operation == ">" && val > value) return true;
    if (operation == "<" && val < value) return true;
 
    return false;
}
 
vector<string> TRestAnalysisTree::GetCutObservables(const string& cut_str) {
    const std::vector<string> validOperators = {"==", "<=", ">=", "!=", "=", ">", "<"};
 
    std::vector<string> cuts = Split(cut_str, "&&", false, true);
 
    vector<string> obsNames;
    for (auto& cut : cuts) {
        string operation, observable;
        for (const auto& validOperator : validOperators) {
            if (cut.find(validOperator) != string::npos) {
                obsNames.push_back((string)cut.substr(0, cut.find(validOperator)));
                break;
            }
        }
    }
    return obsNames;
}
 
void TRestAnalysisTree::EnableBranches(vector<string> obsNames) {
    for (auto& obsName : obsNames) this->SetBranchStatus(obsName.c_str(), true);
}
 
void TRestAnalysisTree::DisableBranches(vector<string> obsNames) {
    for (auto& obsName : obsNames) this->SetBranchStatus(obsName.c_str(), false);
}
 
void TRestAnalysisTree::EnableAllBranches() { this->SetBranchStatus("*", true); }
 
void TRestAnalysisTree::DisableAllBranches() { this->SetBranchStatus("*", false); }
 
void TRestAnalysisTree::EnableQuickObservableValueSetting() { this->fQuickSetObservableValue = true; }
 
void TRestAnalysisTree::DisableQuickObservableValueSetting() { this->fQuickSetObservableValue = false; }
 
Double_t TRestAnalysisTree::GetObservableIntegral(const TString& obsName, Double_t xLow, Double_t xHigh) {
    Int_t id = GetObservableID((std::string)obsName);
 
    if (id < 0) {
        RESTError << "TRestAnalysisTree::GetObservableIntegral. Observable not found! : " << obsName
                  << RESTendl;
        return 0;
    }
 
    Double_t sum = 0;
    for (Int_t n = 0; n < TTree::GetEntries(); n++) {
        TTree::GetEntry(n);
        Double_t value = GetDblObservableValue(id);
 
        if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue;
        sum += value;
    }
 
    return sum;
}
 
Double_t TRestAnalysisTree::GetObservableAverage(const TString& obsName, Double_t xLow, Double_t xHigh) {
    Int_t id = GetObservableID((std::string)obsName);
 
    if (id < 0) {
        RESTError << "TRestAnalysisTree::GetObservableAverage. Observable not found! : " << obsName
                  << RESTendl;
        return 0;
    }
 
    Double_t sum = 0;
    Int_t N = 0;
    for (Int_t n = 0; n < TTree::GetEntries(); n++) {
        TTree::GetEntry(n);
        Double_t value = GetDblObservableValue(id);
 
        if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue;
        N++;
        sum += value;
    }
 
    if (N <= 0) return 0;
 
    return sum / N;
}
 
Double_t TRestAnalysisTree::GetObservableRMS(const TString& obsName, Double_t xLow, Double_t xHigh) {
    Double_t mean = GetObservableAverage(obsName, xLow, xHigh);
 
    Int_t id = GetObservableID((std::string)obsName);
 
    Double_t sum = 0;
    Int_t N = 0;
    for (Int_t n = 0; n < TTree::GetEntries(); n++) {
        TTree::GetEntry(n);
        Double_t value = GetDblObservableValue(id);
 
        if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue;
        N++;
 
        sum += (value - mean) * (value - mean);
    }
 
    if (N <= 0) return 0;
 
    return TMath::Sqrt(sum / N);
}
 
Double_t TRestAnalysisTree::GetObservableMaximum(const TString& obsName, Double_t xLow, Double_t xHigh) {
    Int_t id = GetObservableID((std::string)obsName);
 
    if (id < 0) {
        RESTError << "TRestAnalysisTree::GetObservableMaximum. Observable not found! : " << obsName
                  << RESTendl;
        return 0;
    }
 
    Double_t max = -DBL_MAX;
    for (Int_t n = 0; n < TTree::GetEntries(); n++) {
        TTree::GetEntry(n);
        Double_t value = GetDblObservableValue(id);
 
        if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue;
        if (value > max) max = value;
    }
 
    return max;
}
 
Double_t TRestAnalysisTree::GetObservableMinimum(const TString& obsName, Double_t xLow, Double_t xHigh) {
    Int_t id = GetObservableID((std::string)obsName);
 
    if (id < 0) {
        RESTError << "TRestAnalysisTree::GetObservableMinimum. Observable not found! : " << obsName
                  << RESTendl;
        return 0;
    }
 
    Double_t min = DBL_MAX;
    for (Int_t n = 0; n < TTree::GetEntries(); n++) {
        TTree::GetEntry(n);
        Double_t value = GetDblObservableValue(id);
 
        if (xLow != -1 && xHigh != -1 && (value < xLow || value > xHigh)) continue;
        if (value < min) min = value;
    }
 
    return min;
}
 
Double_t TRestAnalysisTree::GetObservableContour(const TString& obsName, const TString& obsWeight,
                                                 Double_t level, Int_t nBins, Double_t xLow, Double_t xHigh) {
    if (level > 1 || level < 0) {
        RESTWarning << "Level is : " << level << RESTendl;
        RESTWarning << "Level must be between 0 and 1" << RESTendl;
        return 0;
    }
 
    Double_t integral = this->GetIntegral(obsWeight, xLow, xHigh);
 
    TString histDefinition = Form("hc(%5d,%lf,%lf)", nBins, xLow, xHigh);
    if (xHigh == -1)
        this->Draw(obsName + ">>hc", obsWeight);
    else
        this->Draw(obsName + ">>" + histDefinition, obsWeight);
 
    TH1F* htemp = (TH1F*)gPad->GetPrimitive("hc");
 
    Double_t sum = 0;
    for (int i = 0; i < htemp->GetNbinsX(); i++) {
        sum += htemp->GetBinContent(i + 1);
 
        if (sum > level * integral) return htemp->GetBinCenter(i + 1);
    }
    return 0;
}
 
 
std::vector<std::string> TRestAnalysisTree::GetObservableNames() {
    std::vector<std::string> names;
 
    auto branches = GetListOfBranches();
    for (int i = 0; i < branches->GetEntries(); i++) {
        names.push_back((string)branches->At(i)->GetName());
    }
 
    return names;
}
 
Int_t TRestAnalysisTree::WriteAsTTree(const char* name, Int_t option, Int_t bufsize) {
    TTree* tree = new TTree();
 
    char* dest = (char*)tree + sizeof(TObject);
    char* from = (char*)this + sizeof(TObject);
 
    constexpr size_t n = sizeof(TTree) - sizeof(TObject);
 
    // copy the data of the created TTree to a temporary location
    char* temp = (char*)malloc(n);
    memcpy(temp, dest, n);
 
    // copy the data of this TRestAnalysisTree to TTree, and call the new TTree to write
    memcpy(dest, from, n);
    int result = tree->Write(name, option, bufsize);
 
    // restore the data of the created TTree, then we can free it!
    memcpy(dest, temp, n);
    free(temp);
    delete tree;
 
    return result;
}
 
Bool_t TRestAnalysisTree::AddChainFile(const string& _file) {
    auto file = std::unique_ptr<TFile>{TFile::Open(_file.c_str(), "update")};
    if (!file->IsOpen()) {
        RESTWarning << "TRestAnalysisTree::AddChainFile(): failed to open file " << _file << RESTendl;
        return false;
    }
 
    if (this->GetRunOrigin() == 0) {
        this->GetEntry(0);
    }
 
    auto tree = (TRestAnalysisTree*)file->Get("AnalysisTree");
    if (tree != nullptr && tree->GetEntry(0)) {
        if (tree->GetEntry(0) > 0) {
            if (tree->GetRunOrigin() == this->GetRunOrigin()) {
                // this is a valid tree
                delete tree;
 
                if (fChain == nullptr) {
                    fChain = new TChain("AnalysisTree", "Chained AnalysisTree");
                    fChain->Add(this->GetCurrentFile()->GetName());
                }
                return fChain->Add(_file.c_str()) == 1;
            }
            RESTWarning
                << "TRestAnalysisTree::AddChainFile(): invalid file, AnalysisTree in file has different "
                   "run id!"
                << RESTendl;
        }
        RESTWarning << "TRestAnalysisTree::AddChainFile(): invalid file, AnalysisTree in file is empty!"
                    << RESTendl;
        delete tree;
    }
    RESTWarning << "TRestAnalysisTree::AddChainFile(): invalid file, file does not contain an AnalysisTree!"
                << RESTendl;
 
    return false;
}
 
TTree* TRestAnalysisTree::GetTree() const {
    if (fChain == nullptr) {
        return TTree::GetTree();
    }
    return (TTree*)fChain;
}
 
Long64_t TRestAnalysisTree::LoadTree(Long64_t entry) {
    if (fChain == nullptr) {
        return TTree::LoadTree(entry);
    } else {
        return fChain->LoadTree(entry);
    }
 
    return -2;
}
 
Long64_t TRestAnalysisTree::GetEntries() const {
    if (fChain == nullptr) {
        return TTree::GetEntries();
    }
    return fChain->GetEntries();
}
 
Long64_t TRestAnalysisTree::GetEntries(const char* sel) {
    if (fChain == nullptr) {
        return TTree::GetEntries(sel);
    }
    return fChain->GetEntries(sel);
}
 
TRestAnalysisTree::~TRestAnalysisTree() { delete fSubEventTag; }