//Renata Kopecna

#include <bu2kstarmumu_loader.hh>

#include <options.hh>
#include <event.hh>
#include <helpers.hh>

#include <spdlog.h>

#include <TChain.h>
#include <TFile.h>


using namespace fcnc;

///calculate decay angles according to theory convention http://arxiv.org/abs/0805.2525
void bu2kstarmumu_loader::theory_angles(const TLorentzVector& muplus, const TLorentzVector& muminus, const TLorentzVector& kaon, const TLorentzVector& pion, double& costhetal, double& costhetak, double& phi)
{
    //set up boost vectors
    TLorentzVector b = muplus + muminus + kaon + pion;
    TLorentzVector mumu = muplus + muminus;
    TLorentzVector kpi = kaon + pion;
    TVector3 mumuboost(-mumu.BoostVector());
    TVector3 kpiboost(-kpi.BoostVector());
    TVector3 bboost(-b.BoostVector());

    //determine costhetal
    TLorentzVector muminusd(muminus);
    muminusd.Boost(mumuboost);
    TLorentzVector bd(b);
    bd.Boost(mumuboost);
    costhetal = cos(muminusd.Vect().Angle(bd.Vect()));

    //determine costhetak
    TLorentzVector kaondd(kaon);
    kaondd.Boost(kpiboost);
    TLorentzVector bdd(b);
    bdd.Boost(kpiboost);
    costhetak = cos(kaondd.Vect().Angle(bdd.Vect()));

    //determine phi
    TLorentzVector kaonddd(kaon);
    kaonddd.Boost(bboost);
    TLorentzVector pionddd(pion);
    pionddd.Boost(bboost);
    TLorentzVector muminusddd(muminus);
    muminusddd.Boost(bboost);
    TLorentzVector muplusddd(muplus);
    muplusddd.Boost(bboost);
    TVector3 normalkpi = kaonddd.Vect().Cross(pionddd.Vect());
    TVector3 normalmumu = muminusddd.Vect().Cross(muplusddd.Vect());

    phi = normalkpi.Angle(normalmumu);
    TLorentzVector kpiddd(kpi);
    kpiddd.Boost(bboost);
    if (normalkpi.Cross(normalmumu).Dot(kpiddd.Vect()) < 0.0)  phi = -phi;

}

///calculate decay angles according to lhcb convention http://arxiv.org/abs/1304.6325, bzero is true for B0 and false for B0bar
void bu2kstarmumu_loader::lhcb_angles(bool bplus, const TLorentzVector& muplus, const TLorentzVector& muminus, const TLorentzVector& kaon, const TLorentzVector& pion, double& costhetal, double& costhetak, double& phi){

    //set up boost vectors
    TLorentzVector b = muplus + muminus + kaon + pion;
    TLorentzVector mumu = muplus + muminus;
    TLorentzVector kpi = kaon + pion;
    TVector3 mumuboost(-mumu.BoostVector());
    TVector3 kpiboost(-kpi.BoostVector());
    TVector3 bboost(-b.BoostVector());

    //determine costhetal
    TLorentzVector muminusd(muminus);
    muminusd.Boost(mumuboost);
    TLorentzVector muplusd(muplus);
    muplusd.Boost(mumuboost);
    TLorentzVector bd(b);
    bd.Boost(mumuboost);
    if (bplus)    costhetal = cos(muplusd.Vect().Angle(-bd.Vect()));
    else          costhetal = cos(muminusd.Vect().Angle(-bd.Vect()));

    //determine costhetak
    TLorentzVector kaondd(kaon);
    kaondd.Boost(kpiboost);
    TLorentzVector bdd(b);
    bdd.Boost(kpiboost);
    costhetak = cos(kaondd.Vect().Angle(-bdd.Vect()));

    //determine phi
    TLorentzVector kaonddd(kaon);
    kaonddd.Boost(bboost);
    TLorentzVector pionddd(pion);
    pionddd.Boost(bboost);
    TLorentzVector muminusddd(muminus);
    muminusddd.Boost(bboost);
    TLorentzVector muplusddd(muplus);
    muplusddd.Boost(bboost);
    TVector3 normalkpi = kaonddd.Vect().Cross(pionddd.Vect());
    TVector3 normalmumu = muplusddd.Vect().Cross(muminusddd.Vect());

    TLorentzVector kpiddd(kpi);
    kpiddd.Boost(bboost);

    if (bplus)    {
        phi = normalkpi.Angle(normalmumu);
        if ((normalmumu.Cross(normalkpi)).Dot(kpiddd.Vect()) < 0.0)
            phi = -phi;
    }
    else
    {
        phi = normalkpi.Angle(-normalmumu);
        if ((normalmumu.Cross(normalkpi)).Dot(kpiddd.Vect()) < 0.0)
            phi = -phi;
    }

    spdlog::trace("Calculated angles: \t ctl= {0:0.2f},\t ctk = {1:0.2f}, \t phi = {1:0.2f}", costhetal, costhetak ,phi);
    return;
}

//calculate angles (Heidelberg2018). Should be 100% identical to LHCb convention
//definition as in LHCb-ANA-2013-097 page 23
void bu2kstarmumu_loader::hd_angles(bool bplus, const TLorentzVector& muplus, const TLorentzVector& muminus, const TLorentzVector& kaon, const TLorentzVector& pion, double& costhetal_lhs, double& costhetal_rhs, double& costhetak_lhs, double& costhetak_rhs, double& phi){

    //4D and boost vector of B system
    TLorentzVector B = (muplus + muminus + kaon + pion);
    TVector3 boost_B = B.BoostVector();

    //ctk
    TLorentzVector kaon_in_B = kaon;
    kaon_in_B.Boost(-boost_B);
    TLorentzVector pion_in_B = pion;
    pion_in_B.Boost(-boost_B);
    TLorentzVector kstar_in_B = kaon_in_B + pion_in_B;

    TLorentzVector kaon_in_kstar = kaon_in_B;
    kaon_in_kstar.Boost(-kstar_in_B.BoostVector());

    costhetak_lhs = kaon_in_kstar.Vect().Unit().Dot(kstar_in_B.Vect().Unit());

    TLorentzVector B_in_kstar = B;
    B_in_kstar.Boost(-boost_B);
    B_in_kstar.Boost(-kstar_in_B.BoostVector());

    costhetak_rhs = kaon_in_kstar.Vect().Unit().Dot(-B_in_kstar.Vect().Unit());


    //ctl
    TLorentzVector muplus_in_B = bplus ? muplus : muminus;
    muplus_in_B.Boost(-boost_B);
    TLorentzVector muminus_in_B = bplus ? muminus : muplus;
    muminus_in_B.Boost(-boost_B);
    TLorentzVector dimuon_in_B = muplus_in_B + muminus_in_B;

    TLorentzVector muon_in_dimuon = muplus_in_B;
    muon_in_dimuon.Boost(-dimuon_in_B.BoostVector());

    costhetal_lhs = muon_in_dimuon.Vect().Unit().Dot(dimuon_in_B.Vect().Unit());

    TLorentzVector B_in_dimuon = B;
    B_in_dimuon.Boost(-boost_B);
    B_in_dimuon.Boost(-dimuon_in_B.BoostVector());

    costhetal_rhs = muon_in_dimuon.Vect().Unit().Dot(-B_in_dimuon.Vect().Unit());

    //phi
    TVector3 normalkpi =  (kaon_in_B.Vect().Unit().Cross(pion_in_B.Vect().Unit())).Unit();
    TVector3 normalmumu = (muplus_in_B.Vect().Unit().Cross(muminus_in_B.Vect().Unit())).Unit();
    phi = normalkpi.Angle(normalmumu);
    if ((bplus ? 1. : -1.)*(normalmumu.Cross(normalkpi)).Dot(kstar_in_B.Vect().Unit()) < 0.0)
        phi = -phi;
}

///B+ case: no further cuts in the tuple. Read (up to) all events without restrictions
std::vector<event> bu2kstarmumu_loader::read_full_tuple(int year, std::string filename, std::string treename, bool isMC, bool mctruth, bool genLevelMC, int nevents) {

    //Set a year for gen lvl MC by hand
    if(genLevelMC) year = 2017;

    check_year(year);

    //open the root file
    TFile* root_file = new TFile(filename.c_str(), "READ");
    assert(root_file->IsOpen());

    //Read the tree
    TTree* tree = static_cast<TTree*>(root_file->Get(treename.c_str()));
    if (!tree){
        //Check if the tree exists, if not, return an empty vector
        spdlog::warn("Could not get tree {0:s} from file {1:s}", treename, filename);
        std::vector<fcnc::event> blank_vector;
        return blank_vector;
        //assert(0);
    }

    std::vector<event> events;
    event meas;

    Double_t mbu = PDGMASS_B;
    Double_t sigma_mbu = 0.0; //TODO: mbu possibly from the fit?
    Double_t mkstar = PDGMASS_K_STAR;

    Double_t
            muplus_px, muplus_py, muplus_pz, muplus_pe,
            muminus_px, muminus_py, muminus_pz, muminus_pe,
            muplus_px_TRUE, muplus_py_TRUE, muplus_pz_TRUE, muplus_pe_TRUE,
            muminus_px_TRUE, muminus_py_TRUE, muminus_pz_TRUE, muminus_pe_TRUE,
            kaon_px, kaon_py, kaon_pz, kaon_pe,
            pion_px, pion_py, pion_pz, pion_pe,
            kaon_pt, pion_pt, muplus_pt, muminus_pt,
            Ks_pi_plus_px, Ks_pi_plus_py, Ks_pi_plus_pz, Ks_pi_plus_pt, Ks_pi_plus_pe,
            Ks_pi_minus_px, Ks_pi_minus_py, Ks_pi_minus_pz, Ks_pi_minus_pt, Ks_pi_minus_pe,
            gamma1_px, gamma1_py, gamma1_pz, gamma1_pt, gamma1_pe,
            gamma2_px, gamma2_py, gamma2_pz, gamma2_pt, gamma2_pe,
            totalweight, delta_weight;

    Double_t ctl_TT, ctk_TT, phi_TT;

    Int_t b_id = 521, muplus_id = -13, muminus_id = +13;//, kaon_id, pion_id;
    Int_t DDLL;

    UInt_t runno = 0;
    ULong64_t gpstime = 0;
    Short_t polarity = 0;

    if(!genLevelMC){
        tree->SetBranchAddress(opts->DTF ? "B_plus_M_DTF" : "B_plus_MM",&mbu);
        tree->SetBranchAddress(opts->DTF ? "B_plus_MERR_DTF" : "B_plus_MMERR",&sigma_mbu);
        tree->SetBranchAddress(opts->DTF ? "K_star_plus_M_scaled" : "K_star_plus_MM", &mkstar);

        tree->SetBranchAddress("Polarity", &polarity);
        tree->SetBranchAddress("runNumber", &runno);

        tree->SetBranchAddress("B_plus_ThetaL", &ctl_TT);
        tree->SetBranchAddress("B_plus_ThetaK", &ctk_TT);
        tree->SetBranchAddress("B_plus_Phi", &phi_TT);
    }

    if (!mctruth && !genLevelMC){
        tree->SetBranchAddress("GpsTime", &gpstime);

        tree->SetBranchAddress((opts->DTF ? "mu_plus_PX_DTF" : "mu_plus_PX"), &muplus_px);
        tree->SetBranchAddress((opts->DTF ? "mu_plus_PY_DTF" : "mu_plus_PY"), &muplus_py);
        tree->SetBranchAddress((opts->DTF ? "mu_plus_PZ_DTF" : "mu_plus_PZ"), &muplus_pz);
        tree->SetBranchAddress((opts->DTF ? "mu_plus_PE_DTF" : "mu_plus_PE"), &muplus_pe);
        tree->SetBranchAddress((opts->DTF ? "mu_plus_PT_DTF" : "mu_plus_PT"), &muplus_pt);

        tree->SetBranchAddress((opts->DTF ? "mu_minus_PX_DTF" : "mu_minus_PX"), &muminus_px);
        tree->SetBranchAddress((opts->DTF ? "mu_minus_PY_DTF" : "mu_minus_PY"), &muminus_py);
        tree->SetBranchAddress((opts->DTF ? "mu_minus_PZ_DTF" : "mu_minus_PZ"), &muminus_pz);
        tree->SetBranchAddress((opts->DTF ? "mu_minus_PE_DTF" : "mu_minus_PE"), &muminus_pe);
        tree->SetBranchAddress((opts->DTF ? "mu_minus_PT_DTF" : "mu_minus_PT"), &muminus_pt);

        if(opts->KS){  //K_S^0 pi^+ channel

            if(opts->DTF){
                tree->SetBranchAddress("Ks_pi_plus_PX_DTF", &Ks_pi_plus_px);
                tree->SetBranchAddress("Ks_pi_plus_PY_DTF", &Ks_pi_plus_py);
                tree->SetBranchAddress("Ks_pi_plus_PZ_DTF", &Ks_pi_plus_pz);
                tree->SetBranchAddress("Ks_pi_plus_PT_DTF", &Ks_pi_plus_pt);
                tree->SetBranchAddress("Ks_pi_plus_PE_DTF", &Ks_pi_plus_pe);

                tree->SetBranchAddress("Ks_pi_minus_PX_DTF", &Ks_pi_minus_px);
                tree->SetBranchAddress("Ks_pi_minus_PY_DTF", &Ks_pi_minus_py);
                tree->SetBranchAddress("Ks_pi_minus_PZ_DTF", &Ks_pi_minus_pz);
                tree->SetBranchAddress("Ks_pi_minus_PT_DTF", &Ks_pi_minus_pt);
                tree->SetBranchAddress("Ks_pi_minus_PE_DTF", &Ks_pi_minus_pe);
            }
            else{
                tree->SetBranchAddress("K_short_PX", &kaon_px);
                tree->SetBranchAddress("K_short_PY", &kaon_py);
                tree->SetBranchAddress("K_short_PZ", &kaon_pz);
                tree->SetBranchAddress("K_short_PE", &kaon_pe);
                tree->SetBranchAddress("K_short_PT", &kaon_pt);
            }

            tree->SetBranchAddress((opts->DTF ? "pi_plus_PX_DTF" : "pi_plus_PX"), &pion_px);
            tree->SetBranchAddress((opts->DTF ? "pi_plus_PY_DTF" : "pi_plus_PY"), &pion_py);
            tree->SetBranchAddress((opts->DTF ? "pi_plus_PZ_DTF" : "pi_plus_PZ"), &pion_pz);
            tree->SetBranchAddress((opts->DTF ? "pi_plus_PE_DTF" : "pi_plus_PE"), &pion_pe);
            tree->SetBranchAddress((opts->DTF ? "pi_plus_PT_DTF" : "pi_plus_PT"), &pion_pt);

        }
        else{   //K^+pi^0 channel
            tree->SetBranchAddress((opts->DTF ? "K_plus_PX_DTF" : "K_short_PX"), &kaon_px);
            tree->SetBranchAddress((opts->DTF ? "K_plus_PY_DTF" : "K_short_PY"), &kaon_py);
            tree->SetBranchAddress((opts->DTF ? "K_plus_PZ_DTF" : "K_short_PZ"), &kaon_pz);
            tree->SetBranchAddress((opts->DTF ? "K_plus_PE_DTF" : "K_short_PE"), &kaon_pe);
            tree->SetBranchAddress((opts->DTF ? "K_plus_PT_DTF" : "K_short_PT"), &kaon_pt);

            if(opts->DTF){
                tree->SetBranchAddress("gamma1_PX_DTF", &gamma1_px);
                tree->SetBranchAddress("gamma1_PY_DTF", &gamma1_py);
                tree->SetBranchAddress("gamma1_PZ_DTF", &gamma1_pz);
                tree->SetBranchAddress("gamma1_PE_DTF", &gamma1_pe);
                tree->SetBranchAddress("gamma1_PT_DTF", &gamma1_pt);

                tree->SetBranchAddress("gamma2_PX_DTF", &gamma2_px);
                tree->SetBranchAddress("gamma2_PY_DTF", &gamma2_py);
                tree->SetBranchAddress("gamma2_PZ_DTF", &gamma2_pz);
                tree->SetBranchAddress("gamma2_PE_DTF", &gamma2_pe);
                tree->SetBranchAddress("gamma2_PT_DTF", &gamma2_pt);
            }
            else{
                tree->SetBranchAddress("pi_zero_resolved_PX", &pion_px);
                tree->SetBranchAddress("pi_zero_resolved_PY", &pion_py);
                tree->SetBranchAddress("pi_zero_resolved_PZ", &pion_pz);
                tree->SetBranchAddress("pi_zero_resolved_PE", &pion_pe);
                tree->SetBranchAddress("pi_zero_resolved_PT", &pion_pt);
            }
        }
    }
    else if(genLevelMC){
        if(opts->KS){
            tree->SetBranchAddress("KS0_TRUEP_X", &kaon_px);
            tree->SetBranchAddress("KS0_TRUEP_Y", &kaon_py);
            tree->SetBranchAddress("KS0_TRUEP_Z", &kaon_pz);
            tree->SetBranchAddress("KS0_TRUEP_E", &kaon_pe);

            tree->SetBranchAddress("piplus0_TRUEP_X", &pion_px);
            tree->SetBranchAddress("piplus0_TRUEP_Y", &pion_py);
            tree->SetBranchAddress("piplus0_TRUEP_Z", &pion_pz);
            tree->SetBranchAddress("piplus0_TRUEP_E", &pion_pe);

            tree->SetBranchAddress("piplus_TRUEP_X", &Ks_pi_plus_px);
            tree->SetBranchAddress("piplus_TRUEP_Y", &Ks_pi_plus_py);
            tree->SetBranchAddress("piplus_TRUEP_Z", &Ks_pi_plus_pz);
            tree->SetBranchAddress("piplus_TRUEP_E", &Ks_pi_plus_pe);

            tree->SetBranchAddress("piminus_TRUEP_X", &Ks_pi_minus_px);
            tree->SetBranchAddress("piminus_TRUEP_Y", &Ks_pi_minus_py);
            tree->SetBranchAddress("piminus_TRUEP_Z", &Ks_pi_minus_pz);
            tree->SetBranchAddress("piminus_TRUEP_E", &Ks_pi_minus_pe);

        }
        else{
            tree->SetBranchAddress("Kplus_TRUEP_X", &kaon_px);
            tree->SetBranchAddress("Kplus_TRUEP_Y", &kaon_py);
            tree->SetBranchAddress("Kplus_TRUEP_Z", &kaon_pz);
            tree->SetBranchAddress("Kplus_TRUEP_E", &kaon_pe);

            tree->SetBranchAddress("pi0_TRUEP_X", &pion_px);
            tree->SetBranchAddress("pi0_TRUEP_Y", &pion_py);
            tree->SetBranchAddress("pi0_TRUEP_Z", &pion_pz);
            tree->SetBranchAddress("pi0_TRUEP_E", &pion_pe);

            tree->SetBranchAddress("muplus_TRUEP_X", &muplus_px_TRUE);
            tree->SetBranchAddress("muplus_TRUEP_Y", &muplus_py_TRUE);
            tree->SetBranchAddress("muplus_TRUEP_Z", &muplus_pz_TRUE);
            tree->SetBranchAddress("muplus_TRUEP_E", &muplus_pe_TRUE);

            tree->SetBranchAddress("muminus_TRUEP_X", &muminus_px_TRUE);
            tree->SetBranchAddress("muminus_TRUEP_Y", &muminus_py_TRUE);
            tree->SetBranchAddress("muminus_TRUEP_Z", &muminus_pz_TRUE);
            tree->SetBranchAddress("muminus_TRUEP_E", &muminus_pe_TRUE);
        }
    }
    else{ //use MC true information
        if(opts->KS){
            tree->SetBranchAddress("K_short_TRUEP_X", &kaon_px);
            tree->SetBranchAddress("K_short_TRUEP_Y", &kaon_py);
            tree->SetBranchAddress("K_short_TRUEP_Z", &kaon_pz);
            tree->SetBranchAddress("K_short_TRUEP_E", &kaon_pe);

            tree->SetBranchAddress("pi_plus_TRUEP_X", &pion_px);
            tree->SetBranchAddress("pi_plus_TRUEP_Y", &pion_py);
            tree->SetBranchAddress("pi_plus_TRUEP_Z", &pion_pz);
            tree->SetBranchAddress("pi_plus_TRUEP_E", &pion_pe);
        }
        else{
            tree->SetBranchAddress("K_plus_TRUEP_X", &kaon_px);
            tree->SetBranchAddress("K_plus_TRUEP_Y", &kaon_py);
            tree->SetBranchAddress("K_plus_TRUEP_Z", &kaon_pz);
            tree->SetBranchAddress("K_plus_TRUEP_E", &kaon_pe);

            tree->SetBranchAddress("pi_zero_resolved_TRUEP_X", &pion_px);
            tree->SetBranchAddress("pi_zero_resolved_TRUEP_Y", &pion_py);
            tree->SetBranchAddress("pi_zero_resolved_TRUEP_Z", &pion_pz);
            tree->SetBranchAddress("pi_zero_resolved_TRUEP_E", &pion_pe);
        }
    }

    //For any MC, always link the TRUE muon as well:
    if(isMC || mctruth){
        tree->SetBranchAddress("mu_plus_TRUEP_X", &muplus_px_TRUE);
        tree->SetBranchAddress("mu_plus_TRUEP_Y", &muplus_py_TRUE);
        tree->SetBranchAddress("mu_plus_TRUEP_Z", &muplus_pz_TRUE);
        tree->SetBranchAddress("mu_plus_TRUEP_E", &muplus_pe_TRUE);

        tree->SetBranchAddress("mu_minus_TRUEP_X", &muminus_px_TRUE);
        tree->SetBranchAddress("mu_minus_TRUEP_Y", &muminus_py_TRUE);
        tree->SetBranchAddress("mu_minus_TRUEP_Z", &muminus_pz_TRUE);
        tree->SetBranchAddress("mu_minus_TRUEP_E", &muminus_pe_TRUE);
    }

    //particle IDs:
    if(genLevelMC){
        tree->SetBranchAddress("Kst_892_plus_MC_MOTHER_ID", &b_id);
    }
    else{
        tree->SetBranchAddress("B_plus_ID", &b_id);
        tree->SetBranchAddress("mu_minus_ID", &muminus_id);
        tree->SetBranchAddress("mu_plus_ID", &muplus_id);
    }

    //for K_S^0 channel, get DD or LL track information
    if(opts->KS && !genLevelMC)tree->SetBranchAddress("KshortDecayInVeLo", &DDLL);

    //load either sWeights (data) or MC reweights (signal, reference and phase-space MC)
    if(!genLevelMC){
        if(isMC){
            if (opts->KS){
                tree->SetBranchAddress("weight2D_nLongTracks", &totalweight);
                tree->SetBranchAddress("delta_weight2D_nLongTracks", &delta_weight);
            }
            else{
                tree->SetBranchAddress("weight2D_nLongTracks_TM_rndGamma", &totalweight);
                tree->SetBranchAddress("delta_weight2D_nLongTracks_TM_rndGamma", &delta_weight);
            }
        }
        else tree->SetBranchAddress("N_Bplus_sw", &totalweight);
    }

    //determine number of events to load:
    unsigned int nentries = tree->GetEntries();
    if (nevents > 0){
        unsigned int n = TMath::Abs(nevents);
        if(n < nentries)nentries = n;
    }
    spdlog::info("Reading {0:d} events from file {1:s}", nentries ,filename);

    bool DDtracks = strcmp(opts->DDLL.c_str(), "DD") == 0;

    std::vector<int>vCorruptedEvents;
    vCorruptedEvents.clear();


    //loop over input tree:
    for(unsigned int j=0; j < nentries; j++){
        tree->GetEntry(j);
        if (j % (tree->GetEntries()/10) == 0) std::cout << "\r[" << round(double(j*100)/tree->GetEntries()) << "%]" << std::flush;

        //Is it a B+ or a B-?
        bool b_plus = b_id > 0;

        //for genLevelMonteCarlo: define muonID accordingly to B_id:
        if(genLevelMC){
            muminus_id = b_plus ? +13 : -13;
            muplus_id  = b_plus ? -13 : +13;
        }

        //skip events in K_S^0 sample that have incorrect tracktype
        if(opts->KS && !genLevelMC){
            if( DDtracks &&  DDLL)continue;
            if(!DDtracks && !DDLL)continue;
        }

        if(genLevelMC || mctruth){ //MC TRUE samples and GenLevel samples use the MC true information:
            muminus_px = muminus_px_TRUE;
            muminus_py = muminus_py_TRUE;
            muminus_pz = muminus_pz_TRUE;
            muminus_pe = muminus_pe_TRUE;
            muplus_px = muplus_px_TRUE;
            muplus_py = muplus_py_TRUE;
            muplus_pz = muplus_pz_TRUE;
            muplus_pe = muplus_pe_TRUE;
        }

        //initialize 4D vectors for angel calculations
        TLorentzVector muplus, muminus, kaon, pion, muplus_TRUE, muminus_TRUE;

        if(muminus_id == muplus_id){
            spdlog::critical("Cannot have mu+ and mu- to have the same ID: {0:d}", muplus_id);
            spdlog::critical("Event ID = {0:d}");
            assert(0);
        }

        if (muminus_id == +13){
            muminus = TLorentzVector(muminus_px, muminus_py, muminus_pz, muminus_pe);
            muminus_TRUE = TLorentzVector(muminus_px_TRUE, muminus_py_TRUE, muminus_pz_TRUE, muminus_pe_TRUE);
        }
        else if(muminus_id == -13){
            muplus = TLorentzVector(muminus_px, muminus_py, muminus_pz, muminus_pe);
            muplus_TRUE = TLorentzVector(muminus_px_TRUE, muminus_py_TRUE, muminus_pz_TRUE, muminus_pe_TRUE);
        }
        else{
           spdlog::critical("Found invalid muon ID in mu-: {0:d}", muminus_id);
            assert(0);
        }
        if (muplus_id == -13){
            muplus = TLorentzVector(muplus_px, muplus_py, muplus_pz, muplus_pe);
            muplus_TRUE = TLorentzVector(muplus_px_TRUE, muplus_py_TRUE, muplus_pz_TRUE, muplus_pe_TRUE);
        }
        else if(muplus_id == +13){
            muminus = TLorentzVector(muplus_px, muplus_py, muplus_pz, muplus_pe);
            muminus_TRUE = TLorentzVector(muplus_px_TRUE, muplus_py_TRUE, muplus_pz_TRUE, muplus_pe_TRUE);
        }
        else{
            spdlog::critical("Found invalid muon ID in mu-: {0:d}", muplus_id);
            assert(0);
        }

        if(opts->DTF && !mctruth && !genLevelMC){
            //If no true information, re-reconstruct the pion+kaon pair
            if(opts->KS){
                pion = TLorentzVector(pion_px, pion_py, pion_pz, pion_pe);
                kaon = TLorentzVector(Ks_pi_plus_px,  Ks_pi_plus_py,  Ks_pi_plus_pz,  Ks_pi_plus_pe)
                        + TLorentzVector(Ks_pi_minus_px, Ks_pi_minus_py, Ks_pi_minus_pz, Ks_pi_minus_pe);
            }
            else{
                pion = TLorentzVector(gamma1_px, gamma1_py, gamma1_pz, gamma1_pe)
                        + TLorentzVector(gamma2_px, gamma2_py, gamma2_pz, gamma2_pe);
                kaon = TLorentzVector(kaon_px, kaon_py, kaon_pz, kaon_pe);
            }
        }
        else{ //in case of true info, read kaon and pion directly
            pion = TLorentzVector(pion_px, pion_py, pion_pz, pion_pe);
            kaon = TLorentzVector(kaon_px, kaon_py, kaon_pz, kaon_pe);
        }


        //In case of genLvl MC recalculate B mass
        if(genLevelMC)  mbu = (kaon+pion+muplus+muminus).M();

        //de-bugging
        if (opts->KS){
            spdlog::trace("KS_piplus_PX:  {0:f}", Ks_pi_plus_px );
            spdlog::trace("KS_piplus_PY:  {0:f}", Ks_pi_plus_py );
            spdlog::trace("KS_piplus_PZ:  {0:f}", Ks_pi_plus_pz );
            spdlog::trace("KS_piplus_PE:  {0:f}", Ks_pi_plus_pe );
            spdlog::trace("KS_piminus_PX: {0:f}", Ks_pi_minus_px);
            spdlog::trace("KS_piminus_PY: {0:f}", Ks_pi_minus_py);
            spdlog::trace("KS_piminus_PZ: {0:f}", Ks_pi_minus_pz);
            spdlog::trace("KS_piminus_PE: {0:f}", Ks_pi_minus_pe);
        }

        //print vectors
        if(spdlog::default_logger_raw()->level()==spdlog::level::trace){
            spdlog::debug(b_plus ? "[B+]" : "[B-]");
            muplus.Print();
            muminus.Print();
            pion.Print();
            kaon.Print();
        }

        //make sure that the TRUE_P's of muons in MC are not identical
        if(isMC && opts->reject_identical_muon_TRUEID){
            if(muplus_TRUE.X() == muminus_TRUE.X() &&
                    muplus_TRUE.Y() == muminus_TRUE.Y() &&
                    muplus_TRUE.Z() == muminus_TRUE.Z() &&
                    muplus_TRUE.T() == muminus_TRUE.T()){
                vCorruptedEvents.push_back(j);
                continue;
            }
        }

        //make sure that mup != mun
        if(muplus.X() == muminus.X() &&
                muplus.Y() == muminus.Y() &&
                muplus.Z() == muminus.Z() &&
                muplus.T() == muminus.T()){
            vCorruptedEvents.push_back(j);
            continue;
        }

        //get pseudo_q2 from opening angle of leptons:
        Double_t pseudo_q2 = muplus.Angle(muminus.Vect());
        // Double_t pseudo_q2 = std::max(muminus.Pt(), muplus.Pt());

        //calculate angles
        Double_t ctl_theory,    ctk_theory,     phi_theory;
        Double_t ctl_lhcb,      ctk_lhcb,       phi_lhcb;
        Double_t ctl_converted, ctk_converted,  phi_converted;
        Double_t ctl_hd_lhs,    ctk_hd_lhs,     phi_hd;
        Double_t ctl_hd_rhs,    ctk_hd_rhs;

        theory_angles(        muplus, muminus, kaon, pion, ctl_theory, ctk_theory, phi_theory);
        lhcb_angles  (b_plus, muplus, muminus, kaon, pion, ctl_lhcb,   ctk_lhcb,   phi_lhcb);
        hd_angles    (b_plus, muplus, muminus, kaon, pion, ctl_hd_lhs, ctl_hd_rhs, ctk_hd_rhs, ctk_hd_lhs, phi_hd);

        //check and convert theory convention to LHCb convention
        ctl_converted = b_plus ? ctl_theory : -ctl_theory;
        ctk_converted = -ctk_theory;
        if (b_plus) phi_converted = phi_theory > 0 ? TMath::Pi() - phi_theory : -TMath::Pi() - phi_theory;
        else        phi_converted = phi_theory;

        //convert tupletoolP2VV angles to cos()
        ctl_TT = cos(ctl_TT);
        ctk_TT = cos(ctk_TT);

        spdlog::trace("Event #{0:d}", j);
        spdlog::trace("B" + std::string(b_plus ? "+" : "-") + "\t\t: Theory \tLHCb \t\tConverted \tHD(lhs) \tHD(rhs)\tTupleToolP2VV");
        spdlog::trace("cos(theta_L)\t:{0:0.4f}\t {0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t", ctl_theory, ctl_lhcb,  ctl_converted, ctl_hd_lhs, ctl_hd_rhs, ctl_TT);
        spdlog::trace("cos(theta_k)\t:{0:0.4f}\t {0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t", ctk_theory, ctk_lhcb , ctk_converted, ctk_hd_lhs, ctk_hd_rhs, ctk_TT);
        spdlog::trace("phit:{0:0.4f}\t {0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t\n",  phi_theory, phi_lhcb, phi_converted, phi_hd, phi_TT);


        //test all angles to be identical:
        assert(fabs(ctl_lhcb - ctl_converted) < 1.0e-8);
        assert(fabs(ctk_lhcb - ctk_converted) < 1.0e-8);
        assert(fabs(phi_lhcb - phi_converted) < 1.0e-8);

        assert(fabs(ctl_lhcb - ctl_hd_lhs) < 1.0e-8);
        assert(fabs(ctl_lhcb - ctl_hd_rhs) < 1.0e-8);
        assert(fabs(ctk_lhcb - ctk_hd_lhs) < 1.0e-8);
        assert(fabs(ctk_lhcb - ctk_hd_rhs) < 1.0e-8);
        assert(fabs(phi_lhcb - phi_hd) < 1.0e-8);


        meas.year         = year;
        meas.m            = mbu;
        meas.sigma_m      = sigma_mbu;
        meas.costhetal    = ctl_lhcb;
        meas.costhetak    = ctk_lhcb;
        meas.phi          = phi_lhcb;
        meas.costhetal_fa = ctl_lhcb;
        meas.costhetak_fa = ctk_lhcb;
        meas.phi_fa       = phi_lhcb;
        meas.q2           = (muplus+muminus).Mag2()/1.0e+6;//in gev2
        meas.p2           = (kaon+pion).Mag2()/1.0e+6;//in gev2
        meas.weight       = totalweight;
        meas.delta_weight = delta_weight;
        meas.mkpi         = mkstar;
        meas.cp_conjugate = b_plus;//this is B+
        meas.max_mu_pt    = TMath::Max(muplus.Pt(), muminus.Pt());
        meas.kaon_pt        = kaon.Pt();
        meas.pseudo_q2    = pseudo_q2;
        meas.magnet       = polarity;


        bool ApplyGenLevelTightCuts = false;
        if(genLevelMC && ApplyGenLevelTightCuts){
            if(muplus.Theta() < 0.005 || muplus.Theta() > 0.4)continue;
            if(muminus.Theta() < 0.005 || muminus.Theta() > 0.4)continue;
            if(pion.Theta() < 0.005 || pion.Theta() > 0.4)continue;

            TLorentzVector piplus(Ks_pi_plus_px, Ks_pi_plus_py, Ks_pi_plus_pz, Ks_pi_plus_pe);
            TLorentzVector piminus(Ks_pi_minus_px, Ks_pi_minus_py, Ks_pi_minus_pz, Ks_pi_minus_pe);
            if(piplus.P() <= 1800)continue;
            if(piminus.P() <= 1800)continue;
            if(piplus.Theta() < 0.005 || piplus.Theta() > 0.4)continue;
            if(piminus.Theta() < 0.005 || piminus.Theta() > 0.4)continue;
        }

        events.push_back(meas);

    }
    if(vCorruptedEvents.size() > 0){
        spdlog::warn("Found {0:d} corrupted events with idendical muons: ",vCorruptedEvents.size());
        for(UInt_t c = 0; c < vCorruptedEvents.size(); c++) spdlog::warn("Events at: {0:d}",vCorruptedEvents.at(c));
    }

    root_file->Close();
    delete root_file;
    spdlog::info("DONE reading {0:d} events from file " + filename,events.size());
    return events;

}


