//offline selection of BDT variables for B+->Kst+mumu
//david gerick


#include "GlobalFunctions.hh"
#include "Utils.hpp"
#include "Paths.hpp"


using namespace std;
using namespace RooFit ;



//////////////////////////////////////////////////////////////////////
/// preSelection()
/// loops over all .root-files created by subjobs in ganga and converts into one tree. trigger, kinematic, PID and
/// some more cuts are applied. Q^2-Binning applied, Kshort becomes extra branch for decay inside or outside of
/// VELO. histograms of inv. mass distributions are created and saved in separate .root-file
///
/// runAllYearsAndMagPol()
/// preSelection is executed for 2011 + 2012 as well as MagUp and MagDown of the decay
///
/// runAllSignalData()
/// all data for this decay is evaluated
///
/// runAllSignalMC()
/// runs all MC files for the signal channel
///
/// runAllRefMC()
/// same as above but for reference channel
///
/// runComplete()
/// every available data for all three (sub)decays is evaluated, data and MC for signal and reference channel
///
/// getTCKs()
/// return the exact events numbers for every MagDown and MagUp, Year and TCK configuration
/// 
/// plotHistos()
/// Selected histograms generated in the preSelection() function are plotted and nicely formatted before 
/// saved as .eps file

int copyDataForTMVA(std::string year = "2011", std::string magnet = "down") {

    string string_input = GetInputFile(year,magnet,true,false,false,false,smallSample);

    TChain* tree = new TChain("DecayTree");
    tree->Add(string_input.c_str());

    string output_path = string_input;
    replace(output_path,".root","_tmp.root");

    if(tree->GetEntries() == 0){
        //If it  complains while running, it's just trying to read also from David's preselected files
        coutERROR("No entries found in TTree from tfile " + string_input + ". Exit program!");
        coutERROR("Exit program!");
        return 0;
    }

    TFile *output = new TFile(output_path.c_str(),"RECREATE");
    tree->SetBranchStatus("*",1);

    TTree * NewTree = tree->CloneTree(0);

    bool TMedBKGCAT  = true;
    bool TMed_noPi0  = true;
    bool TMed     = true;
    bool TM_gammas  = true;
    NewTree->Branch("TMedBKGCAT", &TMedBKGCAT, "TMedBKGCAT/b");
    NewTree->Branch("TMed_noPi0", &TMed_noPi0, "TMed_noPi0/b");
    NewTree->Branch("TMed", &TMed, "TMed/b");
    NewTree->Branch("TM_gammas", &TM_gammas, "TM_gammas/b");

    for (int evt = 0; evt < tree->GetEntries(); evt++){
        tree->GetEntry(evt);
        TMedBKGCAT = true;
        TMed_noPi0 = true;
        TMed = true;
        TM_gammas = true;
        NewTree->Fill();
    }

    output->cd();
    NewTree->Write("",TObject::kWriteDelete);

    output->Close();
    cout << "[INFO]\t\t New tree created." << endl;

    //careful, possibly leaky (wouldnẗ life be so much easier if one could just change ROOT tree on the spot...
    string command = "mv " + string_input + " ";
    replace(string_input,".root","_backup.root");
    command += string_input;
    system(command.c_str());
    replace(string_input,"_backup.root",".root");
    command = "mv " + output_path + " " + string_input;
    system(command.c_str());

    return 1;

}

int copyDataForTMVAAll(int Run = 1){

    for (auto& year : yearsData(Run)){
        if (copyDataForTMVA(year,"down") ==0) return 0;
        if (copyDataForTMVA(year,"up") ==0) return 0;
    }
    return 1;

}


//the function to loop over all events and merge them into one file; apply cut-based selection while merging

int preSelection(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false, bool B0 = false, bool K1 = false, bool Inc = false) {

    TChain* tree = 0;    

    //Reference channel and PHSP are both MC data
    if(ReferenceChannel)MC = true;
    if(PHSP)MC = true;
    if(K1)MC = true;
    if(Inc)MC = true;

    if(ReferenceChannel && PHSP){
    coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
	PHSP = false;
    }

    //Kst2Kpluspi0Resolved case
    if(Kst2Kpluspi0Resolved){
        tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
        if (B0 || K1 || Inc){
            tree->Add(GetInputFileBkg(year,magnet,false,ReferenceChannel,B0, K1, Inc).c_str());
            coutDebug("Adding " + GetInputFileBkg(year,magnet,false,ReferenceChannel,B0, K1, Inc));
        }
        else{
            tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
            coutDebug("Adding " + GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample));
        }
    }

    //Kst2Kpluspi0Merged case
    if(Kst2Kpluspi0Merged){
        tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
        cout << "[DEBUG]\t\t Openning " << GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str() << endl;
    }

    //Kst2Kspiplus case
    if(Kst2Kspiplus){
        tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    int N = tree->GetEntries();
    cout << "Old files contain " << N << " events" <<  endl;

    if(N == 0){
        coutERROR("Tree has no entries. Exit!");
        return 0;
    }

    //Disable all branches
    tree->SetBranchStatus("*",0);

    //activate needed branches
    tree->SetBranchStatus("B_plus_*Decision*",1) ;

    tree->SetBranchStatus("*PID*",1) ;
    tree->SetBranchStatus("*ProbNN*",1) ;
    //tree->SetBranchStatus("*GhostProb",1) ;

    tree->SetBranchStatus("nCandidate",1) ;
    tree->SetBranchStatus("*Tracks",1) ;
    tree->SetBranchStatus("nVeloClusters", 1) ;
    tree->SetBranchStatus("nTTClusters", 1) ;
    tree->SetBranchStatus("nITClusters", 1) ;
    tree->SetBranchStatus("nOTClusters", 1) ;
    tree->SetBranchStatus("nSPDHits", 1) ;
    tree->SetBranchStatus("nPV",1) ;
    tree->SetBranchStatus("eventNumber",1) ;
    tree->SetBranchStatus("runNumber",1) ;
    tree->SetBranchStatus("GpsTime",1) ;
    tree->SetBranchStatus("EventInSequence",1) ;
    tree->SetBranchStatus("totCandidates",1) ;
    tree->SetBranchStatus("Polarity",1) ;
    //keep the TCKinfo for KS (2012early/late) samples:
    if(Kst2Kspiplus)// && year == "2012" && PHSP)
      tree->SetBranchStatus("HLT1TCK", 1);

    //activate cone variables
    if(IsDST(year,MC,ReferenceChannel,PHSP) || Inc || (B0 && ReferenceChannel) || K1 ){
        tree->SetBranchStatus("*_cp*",1);
        tree->SetBranchStatus("*deltaEta*",1);
        tree->SetBranchStatus("*deltaPhi*",1);
        tree->SetBranchStatus("*cmult*",1);
        tree->SetBranchStatus("*pxasy*",1);
        tree->SetBranchStatus("*pyasy*",1);
        tree->SetBranchStatus("*pzasy*",1);
        tree->SetBranchStatus("*ptasy*",1);
        tree->SetBranchStatus("*pasy*",1);
    }

    tree->SetBranchStatus("*M",1) ;
    tree->SetBranchStatus("*MERR",1) ;
    tree->SetBranchStatus("*MM",1) ;
    tree->SetBranchStatus("*MMERR",1) ;
    tree->SetBranchStatus( "*P", 1 );
    tree->SetBranchStatus( "*PX", 1 );
    tree->SetBranchStatus( "*PY", 1);
    tree->SetBranchStatus( "*PZ", 1);
    tree->SetBranchStatus( "*PE", 1);
    tree->SetBranchStatus( "*PT", 1 );
    tree->SetBranchStatus( "*TAU", 1 );
    tree->SetBranchStatus( "*ETA", 1 );   //This should be fixed now
    tree->SetBranchStatus( "*PHI", 1 );
    tree->SetBranchStatus( "*FD*", 1 );
    tree->SetBranchStatus( "*IP*", 1 );

    tree->SetBranchStatus( "*IPCHI2_OWNPV", 1 );
    tree->SetBranchStatus( "*FDCHI2_OWNPV",1 );
    tree->SetBranchStatus( "*OWNPV*", 1 );
    tree->SetBranchStatus( "*DIRA_OWNPV",1);
    tree->SetBranchStatus( "*ENDVERTEX_CHI2",1 );
    tree->SetBranchStatus( "*PV_Z*",1 );
    tree->SetBranchStatus( "*ENDVERTEX*",1 );
    tree->SetBranchStatus( "*DOCA*",1 ); //This should be fixed now
    tree->SetBranchStatus( "*TRACK_CHI2NDOF",1) ;
    tree->SetBranchStatus( "*TRACK_Type",1) ;

    tree->SetBranchStatus( "*chi2",1 );
    tree->SetBranchStatus( "*nDOF",1 );
    tree->SetBranchStatus( "*status",1 );
    tree->SetBranchStatus( "*ctau",1 );

    tree->SetBranchStatus( "*isMuon",1 );
    tree->SetBranchStatus( "*hasRich",1 );
    tree->SetBranchStatus( "*PID*",1 );

    tree->SetBranchStatus( "*_ID",1 );

    tree->SetBranchStatus( "*Phi",1 );
    tree->SetBranchStatus( "*ThetaL",1 );
    tree->SetBranchStatus( "*ThetaK",1 );

    if(MC){
        tree->SetBranchStatus("*TRUE*",1) ;
        tree->SetBranchStatus("*BKGCAT",1) ;
        tree->SetBranchStatus("*MOTHER_ID",1) ;
        tree->SetBranchStatus( "*CosTheta", 1 );
    }

    if (Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
        tree->SetBranchStatus("gamma*ProbNN*",0) ; //Turn off gamma ProbNN
        tree->SetBranchStatus("*_CL",1);
        tree->SetBranchStatus("gamma*CaloHypo_E",1); //possibly use for tests, doesn't hurt to include it
        tree->SetBranchStatus("gamma*CaloHypo_X",1);
        tree->SetBranchStatus("gamma*CaloHypo_Y",1);
        tree->SetBranchStatus("gamma*CaloHypo_Z",1);
    }

    //DTF convergence
    if(UseDTF){
        tree->SetBranchStatus("B_plus_DTF_status",1);
        tree->SetBranchStatus("B_plus_DTF_chi2",1);
    }

    //define variables
    ULong64_t nEventNumber;
    Int_t nSubSetNumber; //split the subset into three subsets for BDT test&training, obtimization, application

    TLorentzVector LorVec_B_plus;
    TLorentzVector LorVec_K_star_plus;
    TLorentzVector LorVec_K_plus;
    TLorentzVector LorVec_pi_zero_merged;
    TLorentzVector LorVec_pi_zero_resolved;
    TLorentzVector LorVec_K_short;
    TLorentzVector LorVec_pi_plus;
    TLorentzVector LorVec_KS_pi_plus;
    TLorentzVector LorVec_KS_pi_minus;
    TLorentzVector LorVec_mu_minus;
    TLorentzVector LorVec_mu_plus;
    TLorentzVector LorVec_gamma1;
    TLorentzVector LorVec_gamma2;
    TLorentzVector LorVec_DiGamma;

    //potentially misIDed particles:
    TLorentzVector LorVec_pi_plus_MisIDed;
    TLorentzVector LorVec_mu_plus_MisIDed;
    TLorentzVector LorVec_B_plus_From_DoubleMisIDed;
    TLorentzVector LorVec_Jpsi_From_DoubleMisIDed;
    TLorentzVector LorVec_K_star_plus_From_DoubleMisIDed;

    //recombined particles:
    TLorentzVector LorVec_DiMuon;
    TLorentzVector LorVec_KaonPion;
    TLorentzVector LorVec_ReconB;
    TLorentzVector LorVec_ReconBPiZero;
    TLorentzVector LorVec_KshortMuMu;
    TLorentzVector LorVec_KplusMuMu;

        ///B+
    Double_t B_plus_M     = 0.;
    Double_t B_plus_MM    = 0.;
    Double_t B_plus_PX    = 0.;
    Double_t B_plus_PY    = 0.;
    Double_t B_plus_PZ    = 0.;
    Double_t B_plus_PE    = 0.;
    Double_t B_plus_PT    = 0.;

    Float_t f_B_plus_M_DTF[100];
    Float_t f_B_plus_MERR_DTF[100];
    Double_t B_plus_M_DTF    = 0.;
    Double_t B_plus_MERR_DTF = 0.;
    Double_t B_plus_P_DTF    = 0.;
    Double_t B_plus_PX_DTF   = 0.;
    Double_t B_plus_PY_DTF   = 0.;
    Double_t B_plus_PZ_DTF   = 0.;
    Double_t B_plus_PE_DTF   = 0.;
    Double_t B_plus_PT_DTF   = 0.;
    Double_t B_plus_ETA_DTF  = 0.;
    Double_t B_plus_PHI_DTF  = 0.;

        ///K*+
    Double_t K_star_plus_M      = 0.;
    Double_t K_star_plus_MM     = 0.;
    Double_t K_star_plus_PX     = 0.;
    Double_t K_star_plus_PY     = 0.;
    Double_t K_star_plus_PZ     = 0.;
    Double_t K_star_plus_PE     = 0.;
    Double_t K_star_plus_PT     = 0.;

    Float_t f_K_star_plus_M_DTF[100];

    Double_t K_star_plus_M_DTF  = 0.;
    Double_t K_star_plus_P_DTF = 0.;
    Double_t K_star_plus_PX_DTF = 0.;
    Double_t K_star_plus_PY_DTF = 0.;
    Double_t K_star_plus_PZ_DTF = 0.;
    Double_t K_star_plus_PE_DTF = 0.;
    Double_t K_star_plus_PT_DTF = 0.;
    Double_t K_star_plus_ETA_DTF= 0.;
    Double_t K_star_plus_PHI_DTF= 0.;

        ///K+
    Double_t K_plus_PX      = 0.;
    Double_t K_plus_PY      = 0.;
    Double_t K_plus_PZ      = 0.;
    Double_t K_plus_PE      = 0.;
    Double_t K_plus_PT      = 0.;

    Float_t f_K_plus_PX_DTF[100];
    Float_t f_K_plus_PY_DTF[100];
    Float_t f_K_plus_PZ_DTF[100];
    Float_t f_K_plus_PE_DTF[100];
    Float_t f_K_plus_ID_DTF[100];

    Double_t K_plus_M_DTF  = 0.; //Just pdg mass
    Double_t K_plus_P_DTF  = 0.;
    Double_t K_plus_PX_DTF = 0.;
    Double_t K_plus_PY_DTF = 0.;
    Double_t K_plus_PZ_DTF = 0.;
    Double_t K_plus_PE_DTF = 0.;
    Double_t K_plus_PT_DTF = 0.;
    Double_t K_plus_ETA_DTF= 0.;
    Double_t K_plus_PHI_DTF= 0.;
    Int_t K_plus_ID_DTF = 0;

    Double_t K_plus_ProbNNk = 0.;


        ///pi0 (merged)
    Double_t pi_zero_merged_M  = 0.;
    Double_t pi_zero_merged_PX = 0.;
    Double_t pi_zero_merged_PY = 0.;
    Double_t pi_zero_merged_PZ = 0.;
    Double_t pi_zero_merged_PE = 0.;
    Double_t pi_zero_merged_PT = 0.;

        ///pi0 (resolved)
    Double_t pi_zero_resolved_M      = 0.;
    Double_t pi_zero_resolved_PX     = 0.;
    Double_t pi_zero_resolved_PY     = 0.;
    Double_t pi_zero_resolved_PZ     = 0.;
    Double_t pi_zero_resolved_PE     = 0.;
    Double_t pi_zero_resolved_PT     = 0.;

    Float_t f_pi_zero_resolved_M_DTF[100];

    Double_t pi_zero_resolved_M_DTF    = 0.;
    Double_t pi_zero_resolved_P_DTF    = 0.;
    Double_t pi_zero_resolved_PX_DTF   = 0.;
    Double_t pi_zero_resolved_PY_DTF   = 0.;
    Double_t pi_zero_resolved_PZ_DTF   = 0.;
    Double_t pi_zero_resolved_PT_DTF   = 0.;
    Double_t pi_zero_resolved_PE_DTF   = 0.;
    Double_t pi_zero_resolved_ETA_DTF  = 0.;
    Double_t pi_zero_resolved_PHI_DTF  = 0.;

        ///Kshort
    Double_t K_short_PX     = 0.;
    Double_t K_short_PY     = 0.;
    Double_t K_short_PZ     = 0.;
    Double_t K_short_PE     = 0.;
    Double_t K_short_PT     = 0.;

    Float_t f_K_short_M_DTF[100];

    Double_t K_short_M_DTF      = 0.;
    Double_t K_short_PT_DTF     = 0.;
    Double_t K_short_ETA_DTF    = 0.;
    Double_t K_short_PHI_DTF    = 0.;

        ///pi+ (KS)
    Double_t Ks_pi_plus_PX     = 0.;
    Double_t Ks_pi_plus_PY     = 0.;
    Double_t Ks_pi_plus_PZ     = 0.;
    Double_t Ks_pi_plus_PE     = 0.;
    Double_t Ks_pi_plus_PT     = 0.;

    Float_t f_Ks_pi_plus_PX_DTF[100];
    Float_t f_Ks_pi_plus_PY_DTF[100];
    Float_t f_Ks_pi_plus_PZ_DTF[100];
    Float_t f_Ks_pi_plus_PE_DTF[100];
    Float_t f_Ks_pi_plus_ID_DTF[100];

    Double_t Ks_pi_plus_PX_DTF = 0.;
    Double_t Ks_pi_plus_PY_DTF = 0.;
    Double_t Ks_pi_plus_PZ_DTF = 0.;
    Double_t Ks_pi_plus_PE_DTF = 0.;
    Double_t Ks_pi_plus_PT_DTF = 0.;
    Double_t Ks_pi_plus_ETA_DTF = 0.;
    Double_t Ks_pi_plus_PHI_DTF = 0.;
    Int_t Ks_pi_plus_ID_DTF = 0;

        ///pi- (KS)
    Double_t Ks_pi_minus_PX     = 0.;
    Double_t Ks_pi_minus_PY     = 0.;
    Double_t Ks_pi_minus_PZ     = 0.;
    Double_t Ks_pi_minus_PE     = 0.;
    Double_t Ks_pi_minus_PT     = 0.;

    Float_t f_Ks_pi_minus_PX_DTF[100];
    Float_t f_Ks_pi_minus_PY_DTF[100];
    Float_t f_Ks_pi_minus_PZ_DTF[100];
    Float_t f_Ks_pi_minus_PE_DTF[100];
    Float_t f_Ks_pi_minus_ID_DTF[100];

    Double_t Ks_pi_minus_PX_DTF  = 0.;
    Double_t Ks_pi_minus_PY_DTF  = 0.;
    Double_t Ks_pi_minus_PZ_DTF  = 0.;
    Double_t Ks_pi_minus_PE_DTF  = 0.;
    Double_t Ks_pi_minus_PT_DTF  = 0.;
    Double_t Ks_pi_minus_ETA_DTF = 0.;
    Double_t Ks_pi_minus_PHI_DTF = 0.;
    Int_t Ks_pi_minus_ID_DTF = 0;

        ///pi+
    Double_t pi_plus_PX     = 0.;
    Double_t pi_plus_PY     = 0.;
    Double_t pi_plus_PZ     = 0.;
    Double_t pi_plus_PE     = 0.;
    Double_t pi_plus_PT     = 0.;

    Float_t f_pi_plus_PX_DTF[100];
    Float_t f_pi_plus_PY_DTF[100];
    Float_t f_pi_plus_PZ_DTF[100];
    Float_t f_pi_plus_PE_DTF[100];
    Float_t f_pi_plus_ID_DTF[100];

    Double_t pi_plus_PX_DTF = 0.;
    Double_t pi_plus_PY_DTF = 0.;
    Double_t pi_plus_PZ_DTF = 0.;
    Double_t pi_plus_PE_DTF = 0.;
    Double_t pi_plus_PT_DTF = 0.;
    Double_t pi_plus_ETA_DTF= 0.;
    Double_t pi_plus_PHI_DTF= 0.;
    Int_t pi_plus_ID_DTF = 0;

        ///mu-
    Double_t mu_minus_PX     = 0.;
    Double_t mu_minus_PY     = 0.;
    Double_t mu_minus_PZ     = 0.;
    Double_t mu_minus_PE     = 0.;
    Double_t mu_minus_PT     = 0.;

    Float_t f_mu_minus_PX_DTF[100];
    Float_t f_mu_minus_PY_DTF[100];
    Float_t f_mu_minus_PZ_DTF[100];
    Float_t f_mu_minus_PE_DTF[100];
    Float_t f_mu_minus_ID_DTF[100];

    Double_t mu_minus_M_DTF = 0.;
    Double_t mu_minus_P_DTF = 0.;
    Double_t mu_minus_PX_DTF = 0.;
    Double_t mu_minus_PY_DTF = 0.;
    Double_t mu_minus_PZ_DTF = 0.;
    Double_t mu_minus_PE_DTF = 0.;
    Double_t mu_minus_PT_DTF = 0.;
    Double_t mu_minus_ETA_DTF= 0.;
    Double_t mu_minus_PHI_DTF= 0.;
    Int_t mu_minus_ID_DTF = 0;

        ///mu+
    Double_t mu_plus_PX     = 0.;
    Double_t mu_plus_PY     = 0.;
    Double_t mu_plus_PZ     = 0.;
    Double_t mu_plus_PE     = 0.;
    Double_t mu_plus_PT     = 0.;

    Float_t f_mu_plus_PX_DTF[100];
    Float_t f_mu_plus_PY_DTF[100];
    Float_t f_mu_plus_PZ_DTF[100];
    Float_t f_mu_plus_PE_DTF[100];
    Float_t f_mu_plus_ID_DTF[100];

    Double_t mu_plus_M_DTF = 0.;
    Double_t mu_plus_P_DTF = 0.;
    Double_t mu_plus_PX_DTF = 0.;
    Double_t mu_plus_PY_DTF = 0.;
    Double_t mu_plus_PZ_DTF = 0.;
    Double_t mu_plus_PE_DTF = 0.;
    Double_t mu_plus_PT_DTF = 0.;
    Double_t mu_plus_ETA_DTF= 0.;
    Double_t mu_plus_PHI_DTF= 0.;
    Int_t mu_plus_ID_DTF = 0;

        ///gamma1
    Double_t gamma1_PX     = 0.;
    Double_t gamma1_PY     = 0.;
    Double_t gamma1_PZ     = 0.;
    Double_t gamma1_PE     = 0.;
    Double_t gamma1_PT     = 0.;

    Float_t f_gamma1_PX_DTF[100];
    Float_t f_gamma1_PY_DTF[100];
    Float_t f_gamma1_PZ_DTF[100];
    Float_t f_gamma1_PE_DTF[100];
    Float_t f_gamma1_ID_DTF[100];

    Double_t gamma1_M_DTF = 0.;
    Double_t gamma1_P_DTF = 0.;
    Double_t gamma1_PX_DTF = 0.;
    Double_t gamma1_PY_DTF = 0.;
    Double_t gamma1_PZ_DTF = 0.;
    Double_t gamma1_PE_DTF = 0.;
    Double_t gamma1_PT_DTF = 0.;
    Double_t gamma1_ETA_DTF= 0.;
    Double_t gamma1_PHI_DTF= 0.;
    Int_t gamma1_ID_DTF = 0;

        ///gamma2
    Double_t gamma2_PX     = 0.;
    Double_t gamma2_PY     = 0.;
    Double_t gamma2_PZ     = 0.;
    Double_t gamma2_PE     = 0.;
    Double_t gamma2_PT     = 0.;

    Float_t f_gamma2_PX_DTF[100];
    Float_t f_gamma2_PY_DTF[100];
    Float_t f_gamma2_PZ_DTF[100];
    Float_t f_gamma2_PE_DTF[100];
    Float_t f_gamma2_ID_DTF[100];

    Double_t gamma2_M_DTF = 0.;
    Double_t gamma2_P_DTF = 0.;
    Double_t gamma2_PX_DTF = 0.;
    Double_t gamma2_PY_DTF = 0.;
    Double_t gamma2_PZ_DTF = 0.;
    Double_t gamma2_PE_DTF = 0.;
    Double_t gamma2_PT_DTF = 0.;
    Double_t gamma2_ETA_DTF= 0.;
    Double_t gamma2_PHI_DTF= 0.;
    Int_t gamma2_ID_DTF = 0;

    //vertex quality:
    Double_t B_plus_DIRA_OWNPV        = 0.;
    Double_t B_plus_FDCHI2_OWNPV      = 0.;
    Double_t B_plus_IPCHI2_OWNPV      = 0.;
    Double_t K_star_plus_IPCHI2_OWNPV = 0.;
    Double_t mu_plus_IPCHI2_OWNPV     = 0.;
    Double_t mu_minus_IPCHI2_OWNPV    = 0.;
    Double_t K_star_plus_FDCHI2_OWNPV = 0.;


    //particle identification
    Double_t K_plus_PIDK  = 0.;
    Double_t pi_plus_PIDK = 0.;

    Double_t mu_plus_ProbNNmu  = 0.;
    Double_t mu_minus_ProbNNmu = 0.;

    //track type variables
    Int_t Ks_pi_plus_TRACK_Type;
    Int_t Ks_pi_minus_TRACK_Type;

    //confidence level
    Double_t gamma1_CL = 0.;
    Double_t gamma2_CL = 0.;

    //Angles
    Double_t ctl = 0.;
    Double_t ctk = 0.;
    Double_t phi = 0.;

    //trigger lines

    //L0
    Bool_t B_plus_L0MuonDecision_TOS = false;

    //HLT1
    Bool_t B_plus_Hlt1TrackAllL0Decision_TOS = false;
    Bool_t B_plus_Hlt1TrackMuonDecision_TOS = false;

    //HLT2
    Bool_t B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS = false;
    Bool_t B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS = false;

    Bool_t B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS = false;
    Bool_t B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS = false;

    Bool_t B_plus_Hlt2Topo2BodyBBDTDecision_TOS = false;
    Bool_t B_plus_Hlt2Topo3BodyBBDTDecision_TOS = false;

//    Bool_t B_plus_Hlt2SingleMuonDecision_TOS = false;

    Bool_t B_plus_Hlt2DiMuonDetachedDecision_TOS = false;

    //isMuon variables
    Bool_t mu_minus_isMuon = false;
    Bool_t mu_plus_isMuon = false;

    //hasRich variables
    Bool_t pi_plus_hasRich = false;
    Bool_t K_plus_hasRich = false;

    //DTF status
    Float_t B_plus_DTF_status[100];
    Float_t B_plus_DTF_chi2[100];

    ///////////////////////////////////////
    /// LINK VARIABLES TO BRANCHES
    ///////////////////////////////////////

    //event numbers

    tree -> SetBranchAddress( "eventNumber" , &nEventNumber );


    //trigger decisions
    tree -> SetBranchAddress( "B_plus_L0MuonDecision_TOS" , &B_plus_L0MuonDecision_TOS );

    if(year == "2011" || year == "2012"){
        tree -> SetBranchAddress( "B_plus_Hlt1TrackAllL0Decision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );
    }
    else if(year == "2015" || year == "2016" || year == "2017" || year == "2018"){
        tree -> SetBranchAddress( "B_plus_Hlt1TrackMVADecision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );

        if(year == "2016" || year == "2017" || year == "2018"){ //found to be FALSE for all 2015 data
            tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS );
            tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS );
        }
    }
    else{
        coutERROR("Incorrect year given for the assosiation of HLT1 lines with branches! Year given: " + year);
        return 0;
    }
    tree -> SetBranchAddress( "B_plus_Hlt1TrackMuonDecision_TOS" , &B_plus_Hlt1TrackMuonDecision_TOS );
//    tree -> SetBranchAddress( "B_plus_Hlt2SingleMuonDecision_TOS" , &B_plus_Hlt2SingleMuonDecision_TOS );
    tree -> SetBranchAddress( "B_plus_Hlt2DiMuonDetachedDecision_TOS" , &B_plus_Hlt2DiMuonDetachedDecision_TOS );

    //particle identification (PID)
    tree -> SetBranchAddress( "mu_minus_isMuon", &mu_minus_isMuon );
    tree -> SetBranchAddress( "mu_plus_isMuon", &mu_plus_isMuon );


    if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
        tree -> SetBranchAddress( "K_plus_hasRich", &K_plus_hasRich );
        if(UsePIDgen && PHSP){ //TODO add '&& MC' and '&& ReferenceChannel' once the PIDgen tuples are available for these sets
            tree -> SetBranchAddress( "K_plus_PIDK_corr", &K_plus_PIDK );
            tree -> SetBranchAddress( "K_plus_ProbNNk", &K_plus_ProbNNk );
            tree -> SetBranchAddress( "mu_plus_ProbNNmu_corr", &mu_plus_ProbNNmu );
            tree -> SetBranchAddress( "mu_minus_ProbNNmu_corr", &mu_minus_ProbNNmu );
        }
        else{
            //tree -> SetBranchAddress( "K_plus_PIDK", &K_plus_PIDK ); replaced by probNN
            tree -> SetBranchAddress( "K_plus_ProbNNk", &K_plus_ProbNNk );
            tree -> SetBranchAddress( "mu_plus_ProbNNmu", &mu_plus_ProbNNmu );
            tree -> SetBranchAddress( "mu_minus_ProbNNmu", &mu_minus_ProbNNmu );
        }
    }

    if(Kst2Kspiplus){
        tree -> SetBranchAddress( "pi_plus_hasRich", &pi_plus_hasRich );
        if(UsePIDgen && PHSP){ //TODO add '&& MC' and '&& ReferenceChannel' once the PIDgen tuples are available for these sets
            tree -> SetBranchAddress( "pi_plus_PIDK_corr", &pi_plus_PIDK );
            tree -> SetBranchAddress( "mu_plus_ProbNNmu_corr", &mu_plus_ProbNNmu );
            tree -> SetBranchAddress( "mu_minus_ProbNNmu_corr", &mu_minus_ProbNNmu );
        }
        else{
            tree -> SetBranchAddress( "pi_plus_PIDK", &pi_plus_PIDK );
            tree -> SetBranchAddress( "mu_plus_ProbNNmu", &mu_plus_ProbNNmu );
            tree -> SetBranchAddress( "mu_minus_ProbNNmu", &mu_minus_ProbNNmu );
        }
    }

    //particle masses and momenta
    tree -> SetBranchAddress( "B_plus_M" , &B_plus_M );
    tree -> SetBranchAddress( "B_plus_MM" , &B_plus_MM );
    tree -> SetBranchAddress( "B_plus_PX" , &B_plus_PX );
    tree -> SetBranchAddress( "B_plus_PY" , &B_plus_PY );
    tree -> SetBranchAddress( "B_plus_PZ" , &B_plus_PZ );
    tree -> SetBranchAddress( "B_plus_PE" , &B_plus_PE );
    tree -> SetBranchAddress( "B_plus_PT" , &B_plus_PT );
    tree -> SetBranchAddress( "B_plus_DTF_M" , &f_B_plus_M_DTF);
    tree -> SetBranchAddress( "B_plus_DTF_MERR", &f_B_plus_MERR_DTF);

    tree -> SetBranchAddress( "K_star_plus_M" , &K_star_plus_M );
    tree -> SetBranchAddress( "K_star_plus_MM" , &K_star_plus_MM );
    tree -> SetBranchAddress( "K_star_plus_PX" , &K_star_plus_PX );
    tree -> SetBranchAddress( "K_star_plus_PY" , &K_star_plus_PY );
    tree -> SetBranchAddress( "K_star_plus_PZ" , &K_star_plus_PZ );
    tree -> SetBranchAddress( "K_star_plus_PE" , &K_star_plus_PE );
    tree -> SetBranchAddress( "K_star_plus_PT" , &K_star_plus_PT );
    tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_M" , &f_K_star_plus_M_DTF );

    if(Kst2Kpluspi0Merged){
        tree -> SetBranchAddress( "pi_zero_merged_PX" , &pi_zero_merged_PX );
        tree -> SetBranchAddress( "pi_zero_merged_PY" , &pi_zero_merged_PY );
        tree -> SetBranchAddress( "pi_zero_merged_PZ" , &pi_zero_merged_PZ );
        tree -> SetBranchAddress( "pi_zero_merged_PE" , &pi_zero_merged_PE );
        tree -> SetBranchAddress( "pi_zero_merged_PT" , &pi_zero_merged_PT );
        tree -> SetBranchAddress( "pi_zero_merged_M"  , &pi_zero_merged_M  );
    }
    if(Kst2Kpluspi0Resolved){
        tree -> SetBranchAddress( "pi_zero_resolved_PX" , &pi_zero_resolved_PX );
        tree -> SetBranchAddress( "pi_zero_resolved_PY" , &pi_zero_resolved_PY );
        tree -> SetBranchAddress( "pi_zero_resolved_PZ" , &pi_zero_resolved_PZ );
        tree -> SetBranchAddress( "pi_zero_resolved_PE" , &pi_zero_resolved_PE );
        tree -> SetBranchAddress( "pi_zero_resolved_PT" , &pi_zero_resolved_PT );
        tree -> SetBranchAddress( "pi_zero_resolved_M"  , &pi_zero_resolved_M  );
        tree -> SetBranchAddress( "B_plus_DTF_pi0_M"  , &f_pi_zero_resolved_M_DTF  );

        tree -> SetBranchAddress( "gamma1_PX" , &gamma1_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PX" , &f_gamma1_PX_DTF );
        tree -> SetBranchAddress( "gamma1_PY" , &gamma1_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PY" , &f_gamma1_PY_DTF );
        tree -> SetBranchAddress( "gamma1_PZ" , &gamma1_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PZ" , &f_gamma1_PZ_DTF );
        tree -> SetBranchAddress( "gamma1_PE" , &gamma1_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PE" , &f_gamma1_PE_DTF );
        tree -> SetBranchAddress( "gamma1_PT" , &gamma1_PT );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_ID" , &f_gamma1_ID_DTF );
        tree -> SetBranchAddress( "gamma2_PX" , &gamma2_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PX" , &f_gamma2_PX_DTF );
        tree -> SetBranchAddress( "gamma2_PY" , &gamma2_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PY" , &f_gamma2_PY_DTF );
        tree -> SetBranchAddress( "gamma2_PZ" , &gamma2_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PZ" , &f_gamma2_PZ_DTF );
        tree -> SetBranchAddress( "gamma2_PE" , &gamma2_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PE" , &f_gamma2_PE_DTF );
        tree -> SetBranchAddress( "gamma2_PT" , &gamma2_PT );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_ID" , &f_gamma2_ID_DTF );

        tree -> SetBranchAddress( "gamma1_CL" , &gamma1_CL );
        tree -> SetBranchAddress( "gamma2_CL" , &gamma2_CL );
    }
    if(Kst2Kpluspi0Merged || Kst2Kpluspi0Resolved){
        tree -> SetBranchAddress( "K_plus_PX" , &K_plus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PX" , &f_K_plus_PX_DTF );
        tree -> SetBranchAddress( "K_plus_PY" , &K_plus_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PY" , &f_K_plus_PY_DTF );
        tree -> SetBranchAddress( "K_plus_PZ" , &K_plus_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PZ" , &f_K_plus_PZ_DTF );
        tree -> SetBranchAddress( "K_plus_PE" , &K_plus_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PE" , &f_K_plus_PE_DTF );
        tree -> SetBranchAddress( "K_plus_PT" , &K_plus_PT );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_ID" , &f_K_plus_ID_DTF );
    }
    if(Kst2Kspiplus){
        tree -> SetBranchAddress( "K_short_PX" , &K_short_PX );
        tree -> SetBranchAddress( "K_short_PY" , &K_short_PY );
        tree -> SetBranchAddress( "K_short_PZ" , &K_short_PZ );
        tree -> SetBranchAddress( "K_short_PE" , &K_short_PE );
        tree -> SetBranchAddress( "K_short_PT" , &K_short_PT );
        tree -> SetBranchAddress( "B_plus_DTF_KS0_M" , &f_K_short_M_DTF );

        tree -> SetBranchAddress( "pi_plus_PX" , &pi_plus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PX" , &f_pi_plus_PX_DTF );
        tree -> SetBranchAddress( "pi_plus_PY" , &pi_plus_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PY" , &f_pi_plus_PY_DTF );
        tree -> SetBranchAddress( "pi_plus_PZ" , &pi_plus_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PZ" , &f_pi_plus_PZ_DTF );
        tree -> SetBranchAddress( "pi_plus_PE" , &pi_plus_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PE" , &f_pi_plus_PE_DTF );
        tree -> SetBranchAddress( "pi_plus_PT" , &pi_plus_PT );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_ID" , &f_pi_plus_ID_DTF );

        tree -> SetBranchAddress( "Ks_pi_minus_PX" , &Ks_pi_minus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PX" , &f_Ks_pi_minus_PX_DTF );
        tree -> SetBranchAddress( "Ks_pi_minus_PY" , &Ks_pi_minus_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PY" , &f_Ks_pi_minus_PY_DTF );
        tree -> SetBranchAddress( "Ks_pi_minus_PZ" , &Ks_pi_minus_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PZ" , &f_Ks_pi_minus_PZ_DTF );
        tree -> SetBranchAddress( "Ks_pi_minus_PE" , &Ks_pi_minus_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PE" , &f_Ks_pi_minus_PE_DTF );
        tree -> SetBranchAddress( "Ks_pi_minus_PT" , &Ks_pi_minus_PT );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_ID" , &f_Ks_pi_minus_ID_DTF );

        tree -> SetBranchAddress( "Ks_pi_plus_PX" , &Ks_pi_plus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PX" , &f_Ks_pi_plus_PX_DTF );
        tree -> SetBranchAddress( "Ks_pi_plus_PY" , &Ks_pi_plus_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PY" , &f_Ks_pi_plus_PY_DTF );
        tree -> SetBranchAddress( "Ks_pi_plus_PZ" , &Ks_pi_plus_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PZ" , &f_Ks_pi_plus_PZ_DTF );
        tree -> SetBranchAddress( "Ks_pi_plus_PE" , &Ks_pi_plus_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PE" , &f_Ks_pi_plus_PE_DTF );
        tree -> SetBranchAddress( "Ks_pi_plus_PT" , &Ks_pi_plus_PT );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_ID" , &f_Ks_pi_plus_ID_DTF );
    }
    tree -> SetBranchAddress( "mu_minus_PX" , &mu_minus_PX );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PX" , &f_mu_minus_PX_DTF );
    tree -> SetBranchAddress( "mu_minus_PY" , &mu_minus_PY );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PY" , &f_mu_minus_PY_DTF );
    tree -> SetBranchAddress( "mu_minus_PZ" , &mu_minus_PZ );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PZ" , &f_mu_minus_PZ_DTF );
    tree -> SetBranchAddress( "mu_minus_PE" , &mu_minus_PE );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PE" , &f_mu_minus_PE_DTF );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_ID" , &f_mu_minus_ID_DTF );
    tree -> SetBranchAddress( "mu_minus_PT" , &mu_minus_PT );

    tree -> SetBranchAddress( "mu_plus_PX" , &mu_plus_PX );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PX" , &f_mu_plus_PX_DTF );
    tree -> SetBranchAddress( "mu_plus_PY" , &mu_plus_PY );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PY" , &f_mu_plus_PY_DTF );
    tree -> SetBranchAddress( "mu_plus_PZ" , &mu_plus_PZ );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PZ" , &f_mu_plus_PZ_DTF );
    tree -> SetBranchAddress( "mu_plus_PE" , &mu_plus_PE );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PE" , &f_mu_plus_PE_DTF );
    tree -> SetBranchAddress( "mu_plus_PT" , &mu_plus_PT );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_ID" , &f_mu_plus_ID_DTF );

    tree -> SetBranchAddress( "B_plus_DTF_status", &B_plus_DTF_status);
    tree -> SetBranchAddress( "B_plus_DTF_chi2", &B_plus_DTF_chi2);

    tree -> SetBranchAddress( "B_plus_DIRA_OWNPV" , &B_plus_DIRA_OWNPV );
    tree -> SetBranchAddress( "B_plus_FDCHI2_OWNPV" , &B_plus_FDCHI2_OWNPV );
    tree -> SetBranchAddress( "B_plus_IPCHI2_OWNPV" , &B_plus_IPCHI2_OWNPV );
    tree -> SetBranchAddress( "K_star_plus_IPCHI2_OWNPV" , &K_star_plus_IPCHI2_OWNPV );
    tree -> SetBranchAddress( "mu_plus_IPCHI2_OWNPV" , &mu_plus_IPCHI2_OWNPV );
    tree -> SetBranchAddress( "mu_minus_IPCHI2_OWNPV" , &mu_minus_IPCHI2_OWNPV );
    tree -> SetBranchAddress( "K_star_plus_FDCHI2_OWNPV" , &K_star_plus_FDCHI2_OWNPV );

    //Angles
    tree -> SetBranchAddress( "B_plus_ThetaL", &ctl);
    tree -> SetBranchAddress( "B_plus_ThetaK", &ctk);
    tree -> SetBranchAddress( "B_plus_Phi", &phi);

    //track types
    if(Kst2Kspiplus){
        tree -> SetBranchAddress( "Ks_pi_plus_TRACK_Type" , &Ks_pi_plus_TRACK_Type );
        tree -> SetBranchAddress( "Ks_pi_minus_TRACK_Type" , &Ks_pi_minus_TRACK_Type );
    }
    //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
    //Create output file

    TFile* output = 0;
    TFile* histos = 0;

    // Kst2Kpluspi0Resolved case
    if(Kst2Kpluspi0Resolved){
        if (B0 || K1 || Inc){
            output = new TFile(GetInputFileBkg(year,magnet,true,ReferenceChannel,B0,K1,Inc).c_str(),"RECREATE");
            histos = new TFile(GetHistosFileBkg(year,magnet,ReferenceChannel,B0,K1,Inc).c_str(),"RECREATE");
        }
        else{
            output = new TFile(GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
            histos = new TFile(GetHistosFile(year,magnet,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
        }
    }
    // Kst2Kpluspi0Merged case
    if(Kst2Kpluspi0Merged){
        output = new TFile(GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
        histos = new TFile(GetHistosFile(year,magnet,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
    }

    // Kst2Kspiplus case
    if(Kst2Kspiplus){
        output = new TFile(GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
        histos = new TFile(GetHistosFile(year,magnet,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
    }

    output->cd();

    TTree * summary_tree = tree->CloneTree(0);

    //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

    int nCounterTrigger = 0;
    int nCounterKinematic = 0;
    int nCounterPrV = 0;
    int nCounterPoV = 0;
    int nCounterPID = 0;
    int nCounterJpsiResonance = 0;
    int nCounterKstarMassWindow = 0;
    int nCounterB0KsmumuMass = 0;
    int nCounterDiMuonMassInQ2Range = 0;
    int n_DoubleMisID = 0;

    histos->cd();
    TH1 * h_JpsiMass = new TH1D("h_DiMuonMass" , "mass spectra of combined muons", 250, 0, 5000);
    h_JpsiMass->GetXaxis()->SetTitle("DiMuonMass [MeV]");
    h_JpsiMass->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_JpsiMassCut = new TH1D("h_DiMuonMassw/oResonances" , "mass spectra of combined muons after resonance cut", 250, 0, 5000);
    h_JpsiMassCut->GetXaxis()->SetTitle("DiMuonMass [MeV]");
    h_JpsiMassCut->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_JpsiMassResonance = new TH1D("h_DiMuonMassOnlyJpsiResonance" , "mass spectra of combined muons inside the Jpsi band", 75, 2500, 4000);
    h_JpsiMassResonance->GetXaxis()->SetTitle("DiMuonMass [MeV]");
    h_JpsiMassResonance->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_Psi2sMassResonance = new TH1D("h_DiMuonMassOnlyPsi2Resonance" , "mass spectra of combined muons inside the Psi2 band", 75, 2500, 4000);
    h_Psi2sMassResonance->GetXaxis()->SetTitle("DiMuonMass [MeV]");
    h_Psi2sMassResonance->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_CharmoniumResonance = new TH1D("h_DiMuonMassCharmoniumResonance" , "mass spectra of combined muons inside the charmonium bands", 75, 2500, 4000);
    h_CharmoniumResonance->GetXaxis()->SetTitle("DiMuonMass [MeV]");
    h_CharmoniumResonance->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_BplusMass_B02KsmumuMass = new TH1D("h_BplusMass_B02KsmumuMass" , "mass spectra of reconstructed B0 from KS and di-muon", 175, 3500, 6000);
    h_BplusMass_B02KsmumuMass->GetXaxis()->SetTitle("Reconstructed B0 mass [MeV]");
    h_BplusMass_B02KsmumuMass->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_BplusMass_DTF = new TH1D("h_BplusMass_DTF" , "Mass spectra from DTF", 125, 4500, 7000);
    h_BplusMass_DTF->GetXaxis()->SetTitle("B^{+} Mass from DTF [MeV]");
    h_BplusMass_DTF->GetYaxis()->SetTitle("Counts / 20MeV");
    TH2 * h_BplusMass_DiMuonMass = new TH2D("h_BplusMass_DiMuonMass" , "2D mass plot after pre-selection", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
    h_BplusMass_DiMuonMass->GetXaxis()->SetTitle("B^{+} mass [MeV]");
    h_BplusMass_DiMuonMass->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
    TH2 * h_BplusMass_DiMuonMassRaw = new TH2D("h_BplusMass_DiMuonMassRaw" , "2D mass plot after stripping", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
    h_BplusMass_DiMuonMassRaw->GetXaxis()->SetTitle("B^{+} mass [MeV]");
    h_BplusMass_DiMuonMassRaw->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
    TH2 * h_BplusMass_DiMuonMassDD = new TH2D("h_BplusMass_DiMuonMassDD" , "2D mass plot after pre-selection (DD tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
    h_BplusMass_DiMuonMassDD->GetXaxis()->SetTitle("B^{+} mass [MeV]");
    h_BplusMass_DiMuonMassDD->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
    TH2 * h_BplusMass_DiMuonMassRawDD = new TH2D("h_BplusMass_DiMuonMassRawDD" , "2D mass plot after stripping (DD tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
    h_BplusMass_DiMuonMassRawDD->GetXaxis()->SetTitle("B^{+} mass [MeV]");
    h_BplusMass_DiMuonMassRawDD->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
    TH2 * h_BplusMass_DiMuonMassLL = new TH2D("h_BplusMass_DiMuonMassLL" , "2D mass plot after pre-selection (LL tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
    h_BplusMass_DiMuonMassLL->GetXaxis()->SetTitle("B^{+} mass [MeV]");
    h_BplusMass_DiMuonMassLL->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
    TH2 * h_BplusMass_DiMuonMassRawLL = new TH2D("h_BplusMass_DiMuonMassRawLL" , "2D mass plot after stripping (LL tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
    h_BplusMass_DiMuonMassRawLL->GetXaxis()->SetTitle("B^{+} mass [MeV]");
    h_BplusMass_DiMuonMassRawLL->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");


    TH1 * h_ReconstructBmassInCharmoniumResonance = new TH1D("h_K*JpsiMassInCharmoniumResonance" , "mass spectra of K* and Jpsi inside the charmonium bands", 75, 5000, 6500);
    h_ReconstructBmassInCharmoniumResonance->GetXaxis()->SetTitle("K*JpsiMass [MeV]");
    h_ReconstructBmassInCharmoniumResonance->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_ReconstructBmassOffResonance = new TH1D("h_K*JpsiMassOutsideOffCharmoniumResonance" , "mass spectra of K* and DiMuon without the charmonium bands", 75, 5000, 6500);
    h_ReconstructBmassOffResonance->GetXaxis()->SetTitle("K*mumuMass [MeV]");
    h_ReconstructBmassOffResonance->GetYaxis()->SetTitle("Counts / 20MeV");

    TH1 * h_KstarMass = new TH1D("h_KstarMass" , "mass spectra of combined kaon and pion", 90, 500, 1400);
    h_KstarMass->GetXaxis()->SetTitle("KstarMass [MeV]");
    h_KstarMass->GetYaxis()->SetTitle("Counts / 10MeV");
    TH1 * h_KstarMasscut = new TH1D("h_KstarMassIn100MeV window" , "mass spectra of combined kaon and pion in 100MeV window", 90, 500, 1400);
    h_KstarMasscut->GetXaxis()->SetTitle("KstarMass [MeV]");
    h_KstarMasscut->GetYaxis()->SetTitle("Counts / 10MeV");
    TH2 * h_KstarMass2D = new TH2D("h_KstarMass2D" , "mass spectrum of Kstarplus", 90, 500, 1400, 90, 500, 1400);
    h_KstarMass2D->GetXaxis()->SetTitle("mass from tree [MeV]");
    h_KstarMass2D->GetYaxis()->SetTitle("mass from K and pi [MeV]");
    TH2 * h_diMuonReconstructedBMass2D = new TH2D("h_diMuonReconstructedBMass2D" , "mass spectrum of dimuon vs B+candidates", 150, 5000, 6500, 150, 2500, 4000);
    TH2 * h_diMuonReconstructedBMass2DFixPiZero = new TH2D("h_diMuonReconstructedBMass2DFixPiZero" , "mass spectrum of dimuon vs B+candidates with fixed m_{#pi^0}", 150, 5000, 6500, 150, 2500, 4000);
    h_diMuonReconstructedBMass2D->GetYaxis()->SetTitle("m_{#mu#mu} [MeV]");
    if (Kst2Kspiplus)
        h_diMuonReconstructedBMass2D->GetXaxis()->SetTitle("m_{K_{s}#pi^{+}#mu#mu} [MeV]");
    else{
        if (Kst2Kpluspi0Merged)
            h_diMuonReconstructedBMass2D->GetXaxis()->SetTitle("m_{K_{+}#pi^{0}#mu#mu} [MeV]");
        else
            h_diMuonReconstructedBMass2D->GetXaxis()->SetTitle("m_{K_{+}#gamma#gamma#mu#mu} [MeV]");
            h_diMuonReconstructedBMass2DFixPiZero->GetXaxis()->SetTitle("m_{K_{+}#pi^{0}#mu#mu} [MeV]");
    }

    //effect of PID cut on angular distribution
    TH1 * h_BeforePIDcuts_ctl = new TH1D("h_BeforePIDcuts_ctl" , "cos(#theta_{L}) before PID cuts; cos(#theta_{L}); entries", 100, -1., +1.);
    TH1 * h_BeforePIDcuts_ctk = new TH1D("h_BeforePIDcuts_ctk" , "cos(#theta_{K}) before PID cuts; cos(#theta_{K}); entries", 100, -1., +1.);
    TH1 * h_BeforePIDcuts_phi = new TH1D("h_BeforePIDcuts_phi" , "#phi before PID cuts; #phi (rad.); entries", 100, -TMath::Pi(), +TMath::Pi());
    TH1 * h_AfterPIDcuts_ctl  = new TH1D("h_AfterPIDcuts_ctl" , "cos(#theta_{L}) after PID cuts; cos(#theta_{L}); entries", 100, -1., +1.);
    TH1 * h_AfterPIDcuts_ctk  = new TH1D("h_AfterPIDcuts_ctk" , "cos(#theta_{K}) after PID cuts; cos(#theta_{K}); entries", 100, -1., +1.);
    TH1 * h_AfterPIDcuts_phi  = new TH1D("h_AfterPIDcuts_phi" , "#phi after PID cuts; #phi (rad.); entries", 100, -TMath::Pi(), +TMath::Pi());

    /// K_s pi+
    TH2 * h_KshortDaughters = new TH2D("h_KshortDaughters" , "track types of the Kshort daughter pions", 10, -0.5, 9.5, 10, -0.5, 9.5);
    h_KshortDaughters->GetXaxis()->SetTitle("pi minus");
    h_KshortDaughters->GetYaxis()->SetTitle("pi plus");

    TH2 * h_BplusM_vs_BplusM_DoubleMisIDed = new TH2D("h_BplusM_vs_BplusM_DoubleMisIDed" , "comparison of B+ mass versus B+ mass (mu+ and pi+ double misIDed", 100, 5000, 6000, 100, 5000, 6000);
    h_BplusM_vs_BplusM_DoubleMisIDed->GetXaxis()->SetTitle("B_plus_M [MeV]");
    h_BplusM_vs_BplusM_DoubleMisIDed->GetYaxis()->SetTitle("B_plus_M (double misIDed) [MeV]");

    TH2 * h_BplusM_vs_JpsiM_DoubleMisIDed = new TH2D("h_BplusM_vs_JpsiM_DoubleMisIDed" , "overlay of J_psi mass and B+ mass (both mu+ and pi+ double misIDed", 50, 5000, 6000, 200, 0, 4000);
    h_BplusM_vs_JpsiM_DoubleMisIDed->GetXaxis()->SetTitle("B_plus_M (double misIDed) [MeV]");
    h_BplusM_vs_JpsiM_DoubleMisIDed->GetYaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");

    TH2 * h_KstarM_vs_JpsiM_DoubleMisIDed = new TH2D("h_KstarM_vs_JpsiM_DoubleMisIDed" , "overlay of J_psi mass and K*+ mass (both mu+ and pi+ double misIDed", 125, 500, 3000, 200, 0, 4000);
    h_KstarM_vs_JpsiM_DoubleMisIDed->GetXaxis()->SetTitle("K_star_mass (double misIDed) [MeV]");
    h_KstarM_vs_JpsiM_DoubleMisIDed->GetYaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");

    TH1 * h_Jpsi_Spectrum_MisIDed_Muon = new TH1D("h_Jpsi_Spectrum_MisIDed_Muon" , "Invariant mass of Jpsi with a pion misIDed as a muon", 200, 0, 4000);
    h_Jpsi_Spectrum_MisIDed_Muon->GetXaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");
    h_Jpsi_Spectrum_MisIDed_Muon->GetYaxis()->SetTitle("Candidates/20MeV");

    TH1 * h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon = new TH1D("h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon" , "Invariant mass of Jpsi with a pion misIDed as a muon (Veto Resonant DiMuons)", 200, 0, 4000);
    h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->GetXaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");
    h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->GetYaxis()->SetTitle("Candidates/20MeV");

    TH1 * h_Q2_Spectrum_MisIDed_Muon = new TH1D("h_Q2_Spectrum_MisIDed_Muon" , "Q2 of non-resonant DiMuons with a pion misIDed as a muon", 100, 0, 20);
    h_Q2_Spectrum_MisIDed_Muon->GetXaxis()->SetTitle("Q2 (double misIDed) [GeV^{2}]");
    h_Q2_Spectrum_MisIDed_Muon->GetYaxis()->SetTitle("Candidates / 0.2 GeV^{2}");


    /// K+ pi0
    TH1 * h_pi0_resolved_candidates_M = new TH1D("h_pi0_resolved_candidates_M" , "#pi^{0} resolved mass spectra", 70, 100, 170);
    h_pi0_resolved_candidates_M->GetXaxis()->SetTitle("#pi^{0} resolved mass [MeV]");
    h_pi0_resolved_candidates_M->GetYaxis()->SetTitle("Candidates/1MeV");

    TH1 * h_pi0_merged_candidates_M = new TH1D("h_pi0_merged_candidates_M" , "#pi^{0} merged mass spectra", 70, 100, 170);
    h_pi0_merged_candidates_M->GetXaxis()->SetTitle("#pi^{0} merged mass [MeV]");
    h_pi0_merged_candidates_M->GetYaxis()->SetTitle("Candidates/1MeV");

    TH1 * h_pi0_merged_candidates_pT = new TH1D("h_pi0_merged_candidates_pT" , "p_{T} spectra of merged pion", 100, 0, 20000);
    h_pi0_merged_candidates_pT->GetXaxis()->SetTitle("p_{T} (#pi^{0} merged) [MeV]");
    h_pi0_merged_candidates_pT->GetYaxis()->SetTitle("Candidates/200MeV");

    TH1 * h_pi0_resolved_candidates_pT = new TH1D("h_pi0_resolved_candidates_pT" , "p_{T} spectra of resolved pion", 100, 0, 20000);
    h_pi0_resolved_candidates_pT->GetXaxis()->SetTitle("p_{T} (#pi^{0} resolved) [MeV]");
    h_pi0_resolved_candidates_pT->GetYaxis()->SetTitle("Candidates/200MeV");

    TH1 * h_KGammaGammaMuMu_mass_before_preselection = new TH1D("h_KGammaGammaMuMu_mass_before_preselection" , "m_{K#gamma#gamma#mu#mu before}", 100, get_cut_B_plus_M_low(year), cut_B_plus_M_high);
    h_KGammaGammaMuMu_mass_before_preselection->GetXaxis()->SetTitle("m_{K#gamma#gamma#mu#mu} [MeV] before");
    h_KGammaGammaMuMu_mass_before_preselection->GetYaxis()->SetTitle(Form("Counts / %.2fMeV", TMath::Abs(cut_B_plus_M_high - get_cut_B_plus_M_low(year))/100));

    TH1 * h_KGammaGammaMuMu_mass_after_preselection= new TH1D("h_KGammaGammaMuMu_mass_after_preselection" , "m_{K#gamma#gamma#mu#mu after}", 100, get_cut_B_plus_M_low(year), cut_B_plus_M_high);
    h_KGammaGammaMuMu_mass_after_preselection->GetXaxis()->SetTitle("m_{K#gamma#gamma#mu#mu} [MeV] after");
    h_KGammaGammaMuMu_mass_after_preselection->GetYaxis()->SetTitle(Form("Counts / %.2fMeV", TMath::Abs(cut_B_plus_M_high - get_cut_B_plus_M_low(year))/100));

    TH1 * h_KPiMuMu_KGammaGammaMuMu_mass= new TH1D("h_KPiMuMu_KGammaGammaMuMu_mass" , "m_{K#pi_{0}#mu#mu}", 100, get_cut_B_plus_M_low(year), cut_B_plus_M_high);
    h_KPiMuMu_KGammaGammaMuMu_mass->GetXaxis()->SetTitle("m_{K#pi_{0}#mu#mu}[MeV]");
    h_KPiMuMu_KGammaGammaMuMu_mass->GetYaxis()->SetTitle(Form("Counts / %.2fMeV", TMath::Abs(cut_B_plus_M_high - get_cut_B_plus_M_low(year))/100));

    TH2 * h_pi0Mass2D = new TH2D("h_pi0Mass2D" , "mass spectrum of #pi^{0} ", 70, 100, 170, 70, 100, 170);
    h_pi0Mass2D->GetXaxis()->SetTitle("mass from tree [MeV]");
    h_pi0Mass2D->GetYaxis()->SetTitle("mass from #gamma #gamma [MeV]");


    //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

    //Define binning in DiMuon mass^2:
    //bin#                          1      2    3   4    5    6    7    8    9     10     11    12    13    14    15    16    17    18
    float f_BinBoundaries[19] = {0.1, 0.98, 1.1, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 11.0, 11.75, 12.5, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0}; //in (GeV/c2)^2
    for(int i = 0; i < 19; i++)f_BinBoundaries[i] *= 1e6; //in (MeV/c2)^2
    //phi  -resonance: bin #2
    //jpsi -resonance: bin #9
    //psi2s-resonance: bin #12
    unsigned int nDiMuonMassBin = 0;
    summary_tree->Branch("nDiMuonMassBin", &nDiMuonMassBin, "nDiMuonMassBin/I");

    Double_t dQ2;
    summary_tree->Branch("Q2", &dQ2, "dQ2/D");

    unsigned int Is_Phi_JPsi_Psi2S = 0;
    summary_tree->Branch("Is_Phi_JPsi_Psi2S", &Is_Phi_JPsi_Psi2S, "Is_Phi_JPsi_Psi2S/I");

    //generate subsets for BDT training, etc. corresponding to the eventnumber
    //modulus(3) of the eventnumber
    summary_tree->Branch("nSubSetNumber", &nSubSetNumber, "nSubSetNumber/I");

    //check, whether the Kshort Decays within the velo
    //this requires both Kshort daughter pions to be tracktype 3
    unsigned int nKshortDecayInVeLo = 0;
    if(Kst2Kspiplus){
        summary_tree->Branch("KshortDecayInVeLo", &nKshortDecayInVeLo, "nKshortDecayInVeLo/I");
    }

    //save invariant masses of certain combinations into TBranches
    Double_t dDiMuonM;
    summary_tree->Branch("DiMuonInvMass", &dDiMuonM, "dDiMuonM/D");
    Double_t dKshortMuMuM;
    if(Kst2Kspiplus)summary_tree->Branch("KshortMuonsInvMass", &dKshortMuMuM, "dKshortMuMuM/D");
    Double_t dKaonPionM;
    summary_tree->Branch("KaonPionInvMass", &dKaonPionM, "dKaonPionM/D");
    Double_t dReconstructedBM;
    summary_tree->Branch("BplusInvMass", &dReconstructedBM, "dReconstructedBM/D");

    Double_t dJpsiMisIDed_M;
    if(Kst2Kspiplus)summary_tree->Branch("JpsiMisIDedM", &dJpsiMisIDed_M, "JpsiMisIDedM/D");

    Double_t dReconstructedBMFixPiZero = 0;
    if (Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){ //Add B+ mass with fixed pi0 mass
        summary_tree->Branch("BplusInvMass_fixPiZero", &dReconstructedBMFixPiZero, "dReconstructedBMFixPiZero/D");
    }

    //Branch with gamma gamma angle
    Double_t open_angle = 0.;
    if (Kst2Kpluspi0Resolved){
        summary_tree->Branch("open_angle", &open_angle, "open_angle/D");
    }

    //get same name branches with pTasymmetry in cones
    //for the muDST of the data and Run 2 MC, the conept is gotton from the RELINFO tupletools
    //while for Run 1 MC (which is full DST) the conept is determined properly by the default tool and created for multiple cone angles
    Double_t B_plus_ConeDeltaETA,           B_plus_ConeP,           B_plus_ConePasym,           B_plus_ConePT,           B_plus_ConePTasym;

    //create new branches for all trees with the 'new' cone variables! this is needed for the TMVA such that they have the same name
    summary_tree->Branch("B_plus_NEW_ConeDeltaETA", &B_plus_ConeDeltaETA, "B_plus_ConeDeltaETA/D");
    summary_tree->Branch("B_plus_NEW_ConeP",        &B_plus_ConeP,        "B_plus_ConeP/D");
    summary_tree->Branch("B_plus_NEW_ConePasym",    &B_plus_ConePasym,    "B_plus_ConePasym/D");
    summary_tree->Branch("B_plus_NEW_ConePT",       &B_plus_ConePT,       "B_plus_ConePT/D");
    summary_tree->Branch("B_plus_NEW_ConePTasym",   &B_plus_ConePTasym,   "B_plus_ConePTasym/D");

    //Rename DTF branches and save them
    summary_tree->Branch("B_plus_M_DTF",        &B_plus_M_DTF,       "B_plus_M_DTF/D");
    summary_tree->Branch("B_plus_MERR_DTF",     &B_plus_MERR_DTF,    "B_plus_MERR_DTF/D");
    summary_tree->Branch("B_plus_P_DTF",        &B_plus_P_DTF,       "B_plus_P_DTF/D");
    summary_tree->Branch("B_plus_PX_DTF",       &B_plus_PX_DTF,      "B_plus_PX_DTF/D");
    summary_tree->Branch("B_plus_PY_DTF",       &B_plus_PY_DTF,      "B_plus_PY_DTF/D");
    summary_tree->Branch("B_plus_PZ_DTF",       &B_plus_PZ_DTF,      "B_plus_PZ_DTF/D");
    summary_tree->Branch("B_plus_PE_DTF",       &B_plus_PE_DTF,      "B_plus_PE_DTF/D");
    summary_tree->Branch("B_plus_PT_DTF",       &B_plus_PT_DTF,      "B_plus_PT_DTF/D");
    summary_tree->Branch("B_plus_ETA_DTF",      &B_plus_ETA_DTF,     "B_plus_ETA_DTF/D");
    summary_tree->Branch("B_plus_PHI_DTF",      &B_plus_PHI_DTF,     "B_plus_PHI_DTF/D");

    summary_tree->Branch("K_star_plus_M_DTF",   &K_star_plus_M_DTF,  "K_star_plus_M_DTF/D");
    summary_tree->Branch("K_star_plus_P_DTF",   &K_star_plus_P_DTF,  "K_star_plus_P_DTF/D");
    summary_tree->Branch("K_star_plus_PX_DTF",  &K_star_plus_PX_DTF, "K_star_plus_PX_DTF/D");
    summary_tree->Branch("K_star_plus_PY_DTF",  &K_star_plus_PY_DTF, "K_star_plus_PY_DTF/D");
    summary_tree->Branch("K_star_plus_PZ_DTF",  &K_star_plus_PZ_DTF, "K_star_plus_PZ_DTF/D");
    summary_tree->Branch("K_star_plus_PE_DTF",  &K_star_plus_PE_DTF, "K_star_plus_PE_DTF/D");
    summary_tree->Branch("K_star_plus_PT_DTF",  &K_star_plus_PT_DTF, "K_star_plus_PT_DTF/D");
    summary_tree->Branch("K_star_plus_ETA_DTF", &K_star_plus_ETA_DTF,"K_star_plus_ETA_DTF/D");
    summary_tree->Branch("K_star_plus_PHI_DTF", &K_star_plus_PHI_DTF,"K_star_plus_PHI_DTF/D");

    summary_tree->Branch("mu_plus_M_DTF",   &mu_plus_M_DTF,   "mu_plus_M_DTF/D");
    summary_tree->Branch("mu_plus_P_DTF",   &mu_plus_P_DTF,   "mu_plus_P_DTF/D");
    summary_tree->Branch("mu_plus_PX_DTF",  &mu_plus_PX_DTF,  "mu_plus_PX_DTF/D");
    summary_tree->Branch("mu_plus_PY_DTF",  &mu_plus_PY_DTF,  "mu_plus_PY_DTF/D");
    summary_tree->Branch("mu_plus_PZ_DTF",  &mu_plus_PZ_DTF,  "mu_plus_PZ_DTF/D");
    summary_tree->Branch("mu_plus_PE_DTF",  &mu_plus_PE_DTF,  "mu_plus_PE_DTF/D");
    summary_tree->Branch("mu_plus_PT_DTF",  &mu_plus_PT_DTF,  "mu_plus_PT_DTF/D");
    summary_tree->Branch("mu_plus_ETA_DTF", &mu_plus_ETA_DTF, "mu_plus_ETA_DTF/D");
    summary_tree->Branch("mu_plus_PHI_DTF", &mu_plus_PHI_DTF, "mu_plus_PHI_DTF/D");
    summary_tree->Branch("mu_plus_ID_DTF",  &mu_plus_ID_DTF,  "mu_plus_ID_DTF/I");

    summary_tree->Branch("mu_minus_M_DTF",   &mu_minus_M_DTF,  "mu_minus_M_DTF/D");
    summary_tree->Branch("mu_minus_P_DTF",   &mu_minus_P_DTF,  "mu_minus_P_DTF/D");
    summary_tree->Branch("mu_minus_PX_DTF",  &mu_minus_PX_DTF,  "mu_minus_PX_DTF/D");
    summary_tree->Branch("mu_minus_PY_DTF",  &mu_minus_PY_DTF,  "mu_minus_PY_DTF/D");
    summary_tree->Branch("mu_minus_PZ_DTF",  &mu_minus_PZ_DTF,  "mu_minus_PZ_DTF/D");
    summary_tree->Branch("mu_minus_PE_DTF",  &mu_minus_PE_DTF,  "mu_minus_PE_DTF/D");
    summary_tree->Branch("mu_minus_PT_DTF",  &mu_minus_PT_DTF,  "mu_minus_PT_DTF/D");
    summary_tree->Branch("mu_minus_ETA_DTF", &mu_minus_ETA_DTF, "mu_minus_ETA_DTF/D");
    summary_tree->Branch("mu_minus_PHI_DTF", &mu_minus_PHI_DTF, "mu_minus_PHI_DTF/D");
    summary_tree->Branch("mu_minus_ID_DTF",  &mu_minus_ID_DTF,  "mu_minus_ID_DTF/I");

    if (Kst2Kpluspi0Resolved){
        summary_tree->Branch("pi_zero_resolved_M_DTF",   &pi_zero_resolved_M_DTF,   "pi_zero_resolved_M_DTF/D");
        summary_tree->Branch("pi_zero_resolved_P_DTF",   &pi_zero_resolved_P_DTF,   "pi_zero_resolved_P_DTF/D");
        summary_tree->Branch("pi_zero_resolved_PX_DTF",  &pi_zero_resolved_PX_DTF,  "pi_zero_resolved_PX_DTF/D");
        summary_tree->Branch("pi_zero_resolved_PY_DTF",  &pi_zero_resolved_PY_DTF,  "pi_zero_resolved_PY_DTF/D");
        summary_tree->Branch("pi_zero_resolved_PZ_DTF",  &pi_zero_resolved_PZ_DTF,  "pi_zero_resolved_PZ_DTF/D");
        summary_tree->Branch("pi_zero_resolved_PE_DTF",  &pi_zero_resolved_PE_DTF,  "pi_zero_resolved_PE_DTF/D");
        summary_tree->Branch("pi_zero_resolved_PT_DTF",  &pi_zero_resolved_PT_DTF,  "pi_zero_resolved_PT_DTF/D");
        summary_tree->Branch("pi_zero_resolved_ETA_DTF", &pi_zero_resolved_ETA_DTF, "pi_zero_resolved_ETA_DTF/D");
        summary_tree->Branch("pi_zero_resolved_PHI_DTF", &pi_zero_resolved_PHI_DTF, "pi_zero_resolved_PHI_DTF/D");

        summary_tree->Branch("K_plus_M_DTF",   &K_plus_M_DTF,   "K_plus_M_DTF/D");
        summary_tree->Branch("K_plus_P_DTF",   &K_plus_P_DTF,   "K_plus_P_DTF/D");
        summary_tree->Branch("K_plus_PX_DTF",  &K_plus_PX_DTF,  "K_plus_PX_DTF/D");
        summary_tree->Branch("K_plus_PY_DTF",  &K_plus_PY_DTF,  "K_plus_PY_DTF/D");
        summary_tree->Branch("K_plus_PZ_DTF",  &K_plus_PZ_DTF,  "K_plus_PZ_DTF/D");
        summary_tree->Branch("K_plus_PE_DTF",  &K_plus_PE_DTF,  "K_plus_PE_DTF/D");
        summary_tree->Branch("K_plus_PT_DTF",  &K_plus_PT_DTF,  "K_plus_PT_DTF/D");
        summary_tree->Branch("K_plus_ETA_DTF", &K_plus_ETA_DTF, "K_plus_ETA_DTF/D");
        summary_tree->Branch("K_plus_PHI_DTF", &K_plus_PHI_DTF, "K_plus_PHI_DTF/D");
        summary_tree->Branch("K_plus_ID_DTF",  &K_plus_ID_DTF,  "K_plus_ID_DTF/I");

        summary_tree->Branch("gamma1_M_DTF",   &gamma1_M_DTF,   "gamma1_M_DTF/D");
        summary_tree->Branch("gamma1_P_DTF",   &gamma1_P_DTF,   "gamma1_P_DTF/D");
        summary_tree->Branch("gamma1_PX_DTF",  &gamma1_PX_DTF,  "gamma1_PX_DTF/D");
        summary_tree->Branch("gamma1_PY_DTF",  &gamma1_PY_DTF,  "gamma1_PY_DTF/D");
        summary_tree->Branch("gamma1_PZ_DTF",  &gamma1_PZ_DTF,  "gamma1_PZ_DTF/D");
        summary_tree->Branch("gamma1_PE_DTF",  &gamma1_PE_DTF,  "gamma1_PE_DTF/D");
        summary_tree->Branch("gamma1_PT_DTF",  &gamma1_PT_DTF,  "gamma1_PT_DTF/D");
        summary_tree->Branch("gamma1_ETA_DTF", &gamma1_ETA_DTF, "gamma1_ETA_DTF/D");
        summary_tree->Branch("gamma1_PHI_DTF", &gamma1_PHI_DTF, "gamma1_PHI_DTF/D");
        summary_tree->Branch("gamma1_ID_DTF",  &gamma1_ID_DTF,  "gamma1_ID_DTF/I");

        summary_tree->Branch("gamma2_M_DTF",   &gamma2_M_DTF,   "gamma2_M_DTF/D");
        summary_tree->Branch("gamma2_P_DTF",   &gamma2_P_DTF,   "gamma2_P_DTF/D");
        summary_tree->Branch("gamma2_PX_DTF",  &gamma2_PX_DTF,  "gamma2_PX_DTF/D");
        summary_tree->Branch("gamma2_PY_DTF",  &gamma2_PY_DTF,  "gamma2_PY_DTF/D");
        summary_tree->Branch("gamma2_PZ_DTF",  &gamma2_PZ_DTF,  "gamma2_PZ_DTF/D");
        summary_tree->Branch("gamma2_PE_DTF",  &gamma2_PE_DTF,  "gamma2_PE_DTF/D");
        summary_tree->Branch("gamma2_PT_DTF",  &gamma2_PT_DTF,  "gamma2_PT_DTF/D");
        summary_tree->Branch("gamma2_ETA_DTF", &gamma2_ETA_DTF, "gamma2_ETA_DTF/D");
        summary_tree->Branch("gamma2_PHI_DTF", &gamma2_PHI_DTF, "gamma2_PHI_DTF/D");
        summary_tree->Branch("gamma2_ID_DTF",  &gamma2_ID_DTF,  "gamma2_ID_DTF/I");
    }

    if(Kst2Kspiplus){
        summary_tree->Branch("K_short_M_DTF",   &K_short_M_DTF,   "K_short_M_DTF/D");
        summary_tree->Branch("K_short_PT_DTF",  &K_short_PT_DTF,  "K_short_PT_DTF/D");
        summary_tree->Branch("K_short_ETA_DTF", &K_short_ETA_DTF, "K_short_ETA_DTF/D");
        summary_tree->Branch("K_short_PHI_DTF", &K_short_PHI_DTF, "K_short_PHI_DTF/D");

        summary_tree->Branch("pi_plus_PX_DTF",  &pi_plus_PX_DTF,  "pi_plus_PX_DTF/D");
        summary_tree->Branch("pi_plus_PY_DTF",  &pi_plus_PY_DTF,  "pi_plus_PY_DTF/D");
        summary_tree->Branch("pi_plus_PZ_DTF",  &pi_plus_PZ_DTF,  "pi_plus_PZ_DTF/D");
        summary_tree->Branch("pi_plus_PT_DTF",  &pi_plus_PT_DTF,  "pi_plus_PT_DTF/D");
        summary_tree->Branch("pi_plus_ETA_DTF", &pi_plus_ETA_DTF, "pi_plus_ETA_DTF/D");
        summary_tree->Branch("pi_plus_PHI_DTF", &pi_plus_PHI_DTF, "pi_plus_PHI_DTF/D");
        summary_tree->Branch("pi_plus_PE_DTF",  &pi_plus_PE_DTF,  "pi_plus_PE_DTF/D");
        summary_tree->Branch("pi_plus_ID_DTF",  &pi_plus_ID_DTF,  "pi_plus_ID_DTF/I");

        summary_tree->Branch("Ks_pi_plus_PX_DTF",  &Ks_pi_plus_PX_DTF,  "Ks_pi_plus_PX_DTF/D");
        summary_tree->Branch("Ks_pi_plus_PY_DTF",  &Ks_pi_plus_PY_DTF,  "Ks_pi_plus_PY_DTF/D");
        summary_tree->Branch("Ks_pi_plus_PZ_DTF",  &Ks_pi_plus_PZ_DTF,  "Ks_pi_plus_PZ_DTF/D");
        summary_tree->Branch("Ks_pi_plus_PT_DTF",  &Ks_pi_plus_PT_DTF,  "Ks_pi_plus_PT_DTF/D");
        summary_tree->Branch("Ks_pi_plus_ETA_DTF", &Ks_pi_plus_ETA_DTF, "Ks_pi_plus_ETA_DTF/D");
        summary_tree->Branch("Ks_pi_plus_PHI_DTF", &Ks_pi_plus_PHI_DTF, "Ks_pi_plus_PHI_DTF/D");
        summary_tree->Branch("Ks_pi_plus_PE_DTF",  &Ks_pi_plus_PE_DTF,  "Ks_pi_plus_PE_DTF/D");
        summary_tree->Branch("Ks_pi_plus_ID_DTF",  &Ks_pi_plus_ID_DTF,  "Ks_pi_plus_ID_DTF/I");

        summary_tree->Branch("Ks_pi_minus_PX_DTF",  &Ks_pi_minus_PX_DTF,  "Ks_pi_minus_PX_DTF/D");
        summary_tree->Branch("Ks_pi_minus_PY_DTF",  &Ks_pi_minus_PY_DTF,  "Ks_pi_minus_PY_DTF/D");
        summary_tree->Branch("Ks_pi_minus_PZ_DTF",  &Ks_pi_minus_PZ_DTF,  "Ks_pi_minus_PZ_DTF/D");
        summary_tree->Branch("Ks_pi_minus_PT_DTF",  &Ks_pi_minus_PT_DTF,  "Ks_pi_minus_PT_DTF/D");
        summary_tree->Branch("Ks_pi_minus_ETA_DTF", &Ks_pi_minus_ETA_DTF, "Ks_pi_minus_ETA_DTF/D");
        summary_tree->Branch("Ks_pi_minus_PHI_DTF", &Ks_pi_minus_PHI_DTF, "Ks_pi_minus_PHI_DTF/D");
        summary_tree->Branch("Ks_pi_minus_PE_DTF",  &Ks_pi_minus_PE_DTF,  "Ks_pi_minus_PE_DTF/D");
        summary_tree->Branch("Ks_pi_minus_ID_DTF",  &Ks_pi_minus_ID_DTF,  "Ks_pi_minus_ID_DTF/I");
    }

    //link variables according to MC or data for all cone variables
    if (IsDST(year,MC,ReferenceChannel, PHSP) || (MC && B0 && ReferenceChannel) ){
        tree -> SetBranchAddress( "B_plus_deltaEta_1.00" , &B_plus_ConeDeltaETA );
        tree -> SetBranchAddress( "B_plus_cp_1.00" ,       &B_plus_ConeP );
        tree -> SetBranchAddress( "B_plus_pasy_1.00" ,     &B_plus_ConePasym );
        tree -> SetBranchAddress( "B_plus_cpt_1.00" ,      &B_plus_ConePT );
        tree -> SetBranchAddress( "B_plus_ptasy_1.00" ,    &B_plus_ConePTasym );
    }
    else{
        tree -> SetBranchAddress( "B_plus_CONEDELTAETA" ,   &B_plus_ConeDeltaETA );
        tree -> SetBranchAddress( "B_plus_CONEP" ,          &B_plus_ConeP );
        tree -> SetBranchAddress( "B_plus_CONEPASYM" ,      &B_plus_ConePasym );
        tree -> SetBranchAddress( "B_plus_CONEPT" ,         &B_plus_ConePT );
        tree -> SetBranchAddress( "B_plus_CONEPTASYM" ,     &B_plus_ConePTasym );
    }

    bool IsBplus = true;

    ///////////////////////////////////////////////////////////////////////////
    ///
    /// loop over events
    ///
    ///////////////////////////////////////////////////////////////////////////

    for(int i=0; i< N; i++){

        //output the different counters!
        if ((0ul == (i % 10000ul) || i + 1 == N) && i != 0) {
            cout << "Read event " << i << "/" << N
                 << ": tTRG true: " << nCounterTrigger << "/" << i  //events passed trigger cut
                 << " In K*+ mass: " << nCounterKstarMassWindow << "/" << i  //events passed K* mass cut
                 << " kinematics: " << nCounterKinematic << "/" << i  //events passed kinematics cut
                 << " PID: " << nCounterPID << "/" << i  //events passed PID cuts
                 << " Q2 Range: " << nCounterDiMuonMassInQ2Range << "/" << i  //events passed Q^2 cut
                 << " not Jpsi: " << nCounterJpsiResonance << "/" << i  //events passed trigger cut
                 << " B02Ksmm: " << nCounterB0KsmumuMass << "/" << nCounterDiMuonMassInQ2Range  //events passed trigger cut
                 << endl;
        }
        tree->GetEntry(i);

        //If use DTF, check for convergence and kill the loop at an early stage
        if (UseDTF && B_plus_DTF_status[0] > cut_DTF_status) continue;
        if (UseDTF && B_plus_DTF_chi2[0] > cut_DTF_chi2)  continue; //Sanity check on chi2
        //use variable to use consistent pion and muon daugthers for the complex conjugated cases (B-)

        //Correlation between B+ and muons/pions:
        //B_plus_ID     mu_minus_ID     mu_plus_ID      KS_piplus_ID        KS_piminus_ID
        //+521          +13             -13             +211                -211
        //-521          -13             +13             -211                +211

        //in DTF only final state particles are assign an ID, therefore, use the charged daughter of K*+ (i.e. K+ or pi+) for identification
        //B+ ID > 0 <==> K+/pi+ ID > 0
        if(Kst2Kspiplus)IsBplus = f_pi_plus_ID_DTF[0] > 0;
        else IsBplus = f_K_plus_ID_DTF[0] > 0;

        //convert Float_t (arrays) to Double_t
        //momenta
        if(Kst2Kspiplus){
            pi_plus_PX_DTF = (Double_t) f_pi_plus_PX_DTF[0];
            pi_plus_PY_DTF = (Double_t) f_pi_plus_PY_DTF[0];
            pi_plus_PZ_DTF = (Double_t) f_pi_plus_PZ_DTF[0];
            pi_plus_PE_DTF = (Double_t) f_pi_plus_PE_DTF[0];
            pi_plus_ID_DTF = std::round(f_pi_plus_ID_DTF[0]);
            if(IsBplus){
                Ks_pi_minus_PX_DTF = (Double_t) f_Ks_pi_plus_PX_DTF[0];
                Ks_pi_minus_PY_DTF = (Double_t) f_Ks_pi_plus_PY_DTF[0];
                Ks_pi_minus_PZ_DTF = (Double_t) f_Ks_pi_plus_PZ_DTF[0];
                Ks_pi_minus_PE_DTF = (Double_t) f_Ks_pi_plus_PE_DTF[0];
                Ks_pi_minus_ID_DTF = std::round(f_Ks_pi_plus_ID_DTF[0]);

                Ks_pi_plus_PX_DTF = (Double_t) f_Ks_pi_minus_PX_DTF[0];
                Ks_pi_plus_PY_DTF = (Double_t) f_Ks_pi_minus_PY_DTF[0];
                Ks_pi_plus_PZ_DTF = (Double_t) f_Ks_pi_minus_PZ_DTF[0];
                Ks_pi_plus_PE_DTF = (Double_t) f_Ks_pi_minus_PE_DTF[0];
                Ks_pi_plus_ID_DTF = std::round(f_Ks_pi_minus_ID_DTF[0]);
            }
            else{
                Ks_pi_plus_PX_DTF = (Double_t) f_Ks_pi_plus_PX_DTF[0];
                Ks_pi_plus_PY_DTF = (Double_t) f_Ks_pi_plus_PY_DTF[0];
                Ks_pi_plus_PZ_DTF = (Double_t) f_Ks_pi_plus_PZ_DTF[0];
                Ks_pi_plus_PE_DTF = (Double_t) f_Ks_pi_plus_PE_DTF[0];
                Ks_pi_plus_ID_DTF = std::round(f_Ks_pi_plus_ID_DTF[0]);

                Ks_pi_minus_PX_DTF = (Double_t) f_Ks_pi_minus_PX_DTF[0];
                Ks_pi_minus_PY_DTF = (Double_t) f_Ks_pi_minus_PY_DTF[0];
                Ks_pi_minus_PZ_DTF = (Double_t) f_Ks_pi_minus_PZ_DTF[0];
                Ks_pi_minus_PE_DTF = (Double_t) f_Ks_pi_minus_PE_DTF[0];
                Ks_pi_minus_ID_DTF = std::round(f_Ks_pi_minus_ID_DTF[0]);
            }
        }
        if(Kst2Kpluspi0Resolved){
            K_plus_PX_DTF = (Double_t) f_K_plus_PX_DTF[0];
            K_plus_PY_DTF = (Double_t) f_K_plus_PY_DTF[0];
            K_plus_PZ_DTF = (Double_t) f_K_plus_PZ_DTF[0];
            K_plus_PE_DTF = (Double_t) f_K_plus_PE_DTF[0];
            K_plus_ID_DTF = std::round(f_K_plus_ID_DTF[0]);

            gamma1_PX_DTF = (Double_t) f_gamma1_PX_DTF[0];
            gamma1_PY_DTF = (Double_t) f_gamma1_PY_DTF[0];
            gamma1_PZ_DTF = (Double_t) f_gamma1_PZ_DTF[0];
            gamma1_PE_DTF = (Double_t) f_gamma1_PE_DTF[0];
            gamma1_ID_DTF = std::round(f_gamma1_ID_DTF[0]);

            gamma2_PX_DTF = (Double_t) f_gamma2_PX_DTF[0];
            gamma2_PY_DTF = (Double_t) f_gamma2_PY_DTF[0];
            gamma2_PZ_DTF = (Double_t) f_gamma2_PZ_DTF[0];
            gamma2_PE_DTF = (Double_t) f_gamma2_PE_DTF[0];
            gamma2_ID_DTF = std::round(f_gamma2_ID_DTF[0]);
        }

        if(IsBplus){
            mu_minus_PX_DTF = (Double_t) f_mu_plus_PX_DTF[0];
            mu_minus_PY_DTF = (Double_t) f_mu_plus_PY_DTF[0];
            mu_minus_PZ_DTF = (Double_t) f_mu_plus_PZ_DTF[0];
            mu_minus_PE_DTF = (Double_t) f_mu_plus_PE_DTF[0];
            mu_minus_ID_DTF = std::round(f_mu_plus_ID_DTF[0]);

            mu_plus_PX_DTF = (Double_t) f_mu_minus_PX_DTF[0];
            mu_plus_PY_DTF = (Double_t) f_mu_minus_PY_DTF[0];
            mu_plus_PZ_DTF = (Double_t) f_mu_minus_PZ_DTF[0];
            mu_plus_PE_DTF = (Double_t) f_mu_minus_PE_DTF[0];
            mu_plus_ID_DTF = std::round(f_mu_minus_ID_DTF[0]);
        }
        else{
            mu_plus_PX_DTF = (Double_t) f_mu_plus_PX_DTF[0];
            mu_plus_PY_DTF = (Double_t) f_mu_plus_PY_DTF[0];
            mu_plus_PZ_DTF = (Double_t) f_mu_plus_PZ_DTF[0];
            mu_plus_PE_DTF = (Double_t) f_mu_plus_PE_DTF[0];
            mu_plus_ID_DTF = std::round(f_mu_plus_ID_DTF[0]);

            mu_minus_PX_DTF = (Double_t) f_mu_minus_PX_DTF[0];
            mu_minus_PY_DTF = (Double_t) f_mu_minus_PY_DTF[0];
            mu_minus_PZ_DTF = (Double_t) f_mu_minus_PZ_DTF[0];
            mu_minus_PE_DTF = (Double_t) f_mu_minus_PE_DTF[0];
            mu_minus_ID_DTF = std::round(f_mu_minus_ID_DTF[0]);
        }

        //convert masses Float_t to Double_t
        B_plus_M_DTF = (Double_t) f_B_plus_M_DTF[0];
        B_plus_MERR_DTF = (Double_t) f_B_plus_MERR_DTF[0];
        K_star_plus_M_DTF = (Double_t) f_K_star_plus_M_DTF[0];
        if(Kst2Kspiplus){
            K_short_M_DTF = (Double_t) f_K_short_M_DTF[0];
        }
        if(Kst2Kpluspi0Resolved){
            pi_zero_resolved_M_DTF = (Double_t) f_pi_zero_resolved_M_DTF[0];
        }

        //define the Lorentz vectors
        if(UseDTF){ //start useDTF for LorentzVectors
            if(UsePDGmIn4Dvectors){
                LorVec_mu_minus.SetXYZM         (mu_minus_PX_DTF,           mu_minus_PY_DTF,            mu_minus_PZ_DTF,            PDGMASS.MU);
                LorVec_mu_plus.SetXYZM          (mu_plus_PX_DTF,            mu_plus_PY_DTF,             mu_plus_PZ_DTF,             PDGMASS.MU);
                if(Kst2Kspiplus){
                    /// Ks, pi+
                    LorVec_pi_plus.SetXYZM          (pi_plus_PX_DTF,            pi_plus_PY_DTF,             pi_plus_PZ_DTF,             PDGMASS.PI_PLUS);
                    LorVec_KS_pi_plus.SetXYZM       (Ks_pi_plus_PX_DTF,         Ks_pi_plus_PY_DTF,          Ks_pi_plus_PZ_DTF,          PDGMASS.PI_PLUS);
                    LorVec_KS_pi_minus.SetXYZM      (Ks_pi_minus_PX_DTF,        Ks_pi_minus_PY_DTF,         Ks_pi_minus_PZ_DTF,         PDGMASS.PI_PLUS);
                    LorVec_pi_plus_MisIDed.SetXYZM  (pi_plus_PX_DTF,            pi_plus_PY_DTF,             pi_plus_PZ_DTF,             PDGMASS.MU);
                    LorVec_mu_plus_MisIDed.SetXYZM  (mu_plus_PX_DTF,            mu_plus_PY_DTF,             mu_plus_PZ_DTF,             PDGMASS.PI_PLUS);
                }
                else{
                    /// K+, pi0
                    LorVec_K_plus.SetXYZM           (K_plus_PX_DTF,             K_plus_PY_DTF,              K_plus_PZ_DTF,              PDGMASS.K_PLUS);
                    LorVec_gamma1.SetXYZM           (gamma1_PX_DTF,             gamma1_PY_DTF,              gamma1_PZ_DTF,              PDGMASS.GAMMA);
                    LorVec_gamma2.SetXYZM           (gamma2_PX_DTF,             gamma2_PY_DTF,              gamma2_PZ_DTF,              PDGMASS.GAMMA);
                }
            }
            else{
                LorVec_mu_minus.SetPxPyPzE        (mu_minus_PX_DTF,           mu_minus_PY_DTF,            mu_minus_PZ_DTF,            mu_minus_PE_DTF);
                LorVec_mu_plus.SetPxPyPzE         (mu_plus_PX_DTF,            mu_plus_PY_DTF,             mu_plus_PZ_DTF,             mu_plus_PE_DTF);
                if(Kst2Kspiplus){
                    /// Ks, pi+
                    LorVec_pi_plus.SetPxPyPzE         (pi_plus_PX_DTF,            pi_plus_PY_DTF,             pi_plus_PZ_DTF,             pi_plus_PE_DTF);
                    LorVec_KS_pi_plus.SetPxPyPzE      (Ks_pi_plus_PX_DTF,         Ks_pi_plus_PY_DTF,          Ks_pi_plus_PZ_DTF,          Ks_pi_plus_PE_DTF);
                    LorVec_KS_pi_minus.SetPxPyPzE     (Ks_pi_minus_PX_DTF,        Ks_pi_minus_PY_DTF,         Ks_pi_minus_PZ_DTF,         Ks_pi_minus_PE_DTF);
                    LorVec_pi_plus_MisIDed.SetXYZM    (pi_plus_PX_DTF,            pi_plus_PY_DTF,             pi_plus_PZ_DTF,             PDGMASS.MU);
                    LorVec_mu_plus_MisIDed.SetXYZM    (mu_plus_PX_DTF,            mu_plus_PY_DTF,             mu_plus_PZ_DTF,             PDGMASS.PI_PLUS);
                }
                else{
                    /// K+, pi0
                    LorVec_K_plus.SetPxPyPzE    (K_plus_PX_DTF,             K_plus_PY_DTF,              K_plus_PZ_DTF,              K_plus_PE_DTF);
                    LorVec_gamma1.SetPxPyPzE    (gamma1_PX_DTF,             gamma1_PY_DTF,              gamma1_PZ_DTF,              gamma1_PE_DTF);
                    LorVec_gamma2.SetPxPyPzE    (gamma2_PX_DTF,             gamma2_PY_DTF,              gamma2_PZ_DTF,              gamma2_PE_DTF);
                }
            }

            //directly recombine the non-stable particles for later use (K*+ not really needed/used for Kst2Kspiplus)
            if(Kst2Kspiplus){
                LorVec_K_short = LorVec_KS_pi_plus + LorVec_KS_pi_minus;
                LorVec_KaonPion = LorVec_K_short + LorVec_pi_plus;
            }
            if(Kst2Kpluspi0Resolved){
                LorVec_DiGamma = LorVec_gamma1 + LorVec_gamma2;
                LorVec_pi_zero_resolved = LorVec_DiGamma;
                LorVec_KaonPion = LorVec_K_plus + LorVec_pi_zero_resolved;
            }
            LorVec_K_star_plus  = LorVec_KaonPion;
            LorVec_DiMuon       = LorVec_mu_minus       + LorVec_mu_plus;
            LorVec_B_plus       = LorVec_K_star_plus    + LorVec_DiMuon;

            //determine transvers momenta PT for all particles:
            B_plus_P_DTF                = LorVec_B_plus.P();
            B_plus_PX_DTF               = LorVec_B_plus.Px();
            B_plus_PY_DTF               = LorVec_B_plus.Py();
            B_plus_PZ_DTF               = LorVec_B_plus.Pz();
            B_plus_PT_DTF               = LorVec_B_plus.Pt();
            B_plus_PE_DTF               = LorVec_B_plus.E();
            B_plus_ETA_DTF              = LorVec_B_plus.Eta();
            B_plus_PHI_DTF              = LorVec_B_plus.Phi();

            K_star_plus_M_DTF           = LorVec_K_star_plus.M();
            K_star_plus_P_DTF           = LorVec_K_star_plus.P();
            K_star_plus_PX_DTF          = LorVec_K_star_plus.Px();
            K_star_plus_PY_DTF          = LorVec_K_star_plus.Py();
            K_star_plus_PZ_DTF          = LorVec_K_star_plus.Pz();
            K_star_plus_PT_DTF          = LorVec_K_star_plus.Pt();
            K_star_plus_PE_DTF          = LorVec_K_star_plus.E();
            K_star_plus_ETA_DTF         = LorVec_K_star_plus.Eta();
            K_star_plus_PHI_DTF         = LorVec_K_star_plus.Phi();

            mu_plus_P_DTF               = LorVec_mu_plus.P();
            mu_plus_M_DTF               = LorVec_mu_plus.M();
            mu_plus_PT_DTF              = LorVec_mu_plus.Pt();
            mu_plus_ETA_DTF             = LorVec_mu_plus.Eta();
            mu_plus_PHI_DTF             = LorVec_mu_plus.Phi();

            mu_minus_P_DTF              = LorVec_mu_minus.P();
            mu_minus_M_DTF              = LorVec_mu_minus.M();
            mu_minus_PT_DTF             = LorVec_mu_minus.Pt();
            mu_minus_ETA_DTF            = LorVec_mu_minus.Eta();
            mu_minus_PHI_DTF            = LorVec_mu_minus.Phi();
            if(Kst2Kspiplus){
                K_short_PT_DTF          = LorVec_K_short.Pt();
                K_short_ETA_DTF         = LorVec_K_short.Eta();
                K_short_PHI_DTF         = LorVec_K_short.Phi();
                Ks_pi_plus_PT_DTF       = LorVec_KS_pi_plus.Pt();
                Ks_pi_plus_ETA_DTF      = LorVec_KS_pi_plus.Eta();
                Ks_pi_plus_PHI_DTF      = LorVec_KS_pi_plus.Phi();
                Ks_pi_minus_PT_DTF      = LorVec_KS_pi_minus.Pt();
                Ks_pi_minus_ETA_DTF     = LorVec_KS_pi_minus.Eta();
                Ks_pi_minus_PHI_DTF     = LorVec_KS_pi_minus.Phi();
                pi_plus_PT_DTF          = LorVec_pi_plus.Pt();
                pi_plus_ETA_DTF         = LorVec_pi_plus.Eta();
                pi_plus_PHI_DTF         = LorVec_pi_plus.Phi();
            }
            if(Kst2Kpluspi0Resolved){
                K_plus_P_DTF            = LorVec_K_plus.P();
                K_plus_M_DTF            = LorVec_K_plus.M();
                K_plus_PT_DTF           = LorVec_K_plus.Pt();
                K_plus_ETA_DTF          = LorVec_K_plus.Eta();
                K_plus_PHI_DTF          = LorVec_K_plus.Phi();

                pi_zero_resolved_P_DTF	= LorVec_pi_zero_resolved.P();
                pi_zero_resolved_PX_DTF	= LorVec_pi_zero_resolved.Px();
                pi_zero_resolved_PY_DTF	= LorVec_pi_zero_resolved.Py();
                pi_zero_resolved_PZ_DTF	= LorVec_pi_zero_resolved.Pz();
                pi_zero_resolved_PT_DTF	= LorVec_pi_zero_resolved.Pt();
                pi_zero_resolved_PE_DTF	= LorVec_pi_zero_resolved.E();
                pi_zero_resolved_ETA_DTF= LorVec_pi_zero_resolved.Eta();
                pi_zero_resolved_PHI_DTF= LorVec_pi_zero_resolved.Phi();

                gamma1_M_DTF            = LorVec_gamma1.M();
                gamma1_P_DTF            = LorVec_gamma1.P();
                gamma1_PT_DTF           = LorVec_gamma1.Pt();
                gamma1_ETA_DTF          = LorVec_gamma1.Eta();
                gamma1_PHI_DTF          = LorVec_gamma1.Phi();

                gamma2_M_DTF            = LorVec_gamma2.M();
                gamma2_P_DTF            = LorVec_gamma2.P();
                gamma2_PT_DTF           = LorVec_gamma2.Pt();
                gamma2_ETA_DTF          = LorVec_gamma2.Eta();
                gamma2_PHI_DTF          = LorVec_gamma2.Phi();
            }
        } //end useDTF for LorentzVectors

        else{ //NON-DTF variables
            if(UsePDGmIn4Dvectors){ //USE PDG masses
                LorVec_B_plus.SetXYZM           (B_plus_PX,             B_plus_PY,              B_plus_PZ,              PDGMASS.B_PLUS);
                LorVec_K_star_plus.SetXYZM      (K_star_plus_PX,        K_star_plus_PY,         K_star_plus_PZ,         PDGMASS.K_STAR_PLUS);
                LorVec_mu_minus.SetXYZM         (mu_minus_PX,           mu_minus_PY,            mu_minus_PZ,            PDGMASS.MU);
                LorVec_mu_plus.SetXYZM          (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             PDGMASS.MU);
                if(Kst2Kspiplus){
                    /// Ks, pi+
                    LorVec_K_short.SetXYZM          (K_short_PX,            K_short_PY,             K_short_PZ,             PDGMASS.K_SHORT);
                    LorVec_pi_plus.SetXYZM          (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             PDGMASS.PI_PLUS);
                    LorVec_KS_pi_plus.SetXYZM       (Ks_pi_plus_PX,         Ks_pi_plus_PY,          Ks_pi_plus_PZ,          PDGMASS.PI_PLUS);
                    LorVec_KS_pi_minus.SetXYZM      (Ks_pi_minus_PX,        Ks_pi_minus_PY,         Ks_pi_minus_PZ,         PDGMASS.PI_PLUS);

                    LorVec_pi_plus_MisIDed.SetXYZM  (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             PDGMASS.MU);
                    LorVec_mu_plus_MisIDed.SetXYZM  (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             PDGMASS.PI_PLUS);
                }
                else{
                    /// K+, pi0
                    LorVec_K_plus.SetXYZM           (K_plus_PX,             K_plus_PY,              K_plus_PZ,              PDGMASS.K_PLUS);
                    LorVec_pi_zero_merged.SetXYZM   (pi_zero_merged_PX,     pi_zero_merged_PY,      pi_zero_merged_PZ,      PDGMASS.PI_ZERO);
                    LorVec_pi_zero_resolved.SetXYZM (pi_zero_resolved_PX,   pi_zero_resolved_PY,    pi_zero_resolved_PZ,    PDGMASS.PI_ZERO);
                    LorVec_gamma1.SetXYZM           (gamma1_PX,             gamma1_PY,              gamma1_PZ,              PDGMASS.GAMMA);
                    LorVec_gamma2.SetXYZM           (gamma2_PX,             gamma2_PY,              gamma2_PZ,              PDGMASS.GAMMA);
                }
            }
            else{
                LorVec_B_plus.SetPxPyPzE        (B_plus_PX,             B_plus_PY,              B_plus_PZ,              B_plus_PE);
                LorVec_K_star_plus.SetPxPyPzE   (K_star_plus_PX,        K_star_plus_PY,         K_star_plus_PZ,         K_star_plus_PE);
                LorVec_mu_minus.SetPxPyPzE      (mu_minus_PX,           mu_minus_PY,            mu_minus_PZ,            mu_minus_PE);
                LorVec_mu_plus.SetPxPyPzE       (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             mu_plus_PE);
                if(Kst2Kspiplus){
                    /// Ks, pi+
                    LorVec_K_short.SetPxPyPzE       (K_short_PX,            K_short_PY,             K_short_PZ,             K_short_PE);
                    LorVec_pi_plus.SetPxPyPzE       (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             pi_plus_PE);
                    LorVec_KS_pi_plus.SetPxPyPzE    (Ks_pi_plus_PX,         Ks_pi_plus_PY,          Ks_pi_plus_PZ,          Ks_pi_plus_PE);
                    LorVec_KS_pi_minus.SetPxPyPzE   (Ks_pi_minus_PX,        Ks_pi_minus_PY,         Ks_pi_minus_PZ,         Ks_pi_minus_PE);

                    LorVec_pi_plus_MisIDed.SetXYZM  (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             PDGMASS.MU);
                    LorVec_mu_plus_MisIDed.SetXYZM  (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             PDGMASS.PI_PLUS);
                }
                else{
                    /// K+, pi0
                    LorVec_K_plus.SetPxPyPzE            (K_plus_PX,             K_plus_PY,              K_plus_PZ,              K_plus_PE);
                    LorVec_pi_zero_merged.SetPxPyPzE	(pi_zero_merged_PX,     pi_zero_merged_PY,      pi_zero_merged_PZ,      pi_zero_merged_PE);
                    LorVec_pi_zero_resolved.SetPxPyPzE  (pi_zero_resolved_PX,   pi_zero_resolved_PY,    pi_zero_resolved_PZ,    pi_zero_resolved_PE);
                    LorVec_gamma1.SetPxPyPzE            (gamma1_PX,             gamma1_PY,              gamma1_PZ,              gamma1_PE);
                    LorVec_gamma2.SetPxPyPzE            (gamma2_PX,             gamma2_PY,              gamma2_PZ,              gamma2_PE);
                }
            }

            //Recombination of 'Parent'-Particles
            LorVec_DiMuon = LorVec_mu_minus + LorVec_mu_plus;
            if(Kst2Kspiplus){
                LorVec_KaonPion = LorVec_K_short + LorVec_pi_plus;
            }
            if(Kst2Kpluspi0Merged){
                LorVec_KaonPion = LorVec_K_plus + LorVec_pi_zero_merged;
            }
            if(Kst2Kpluspi0Resolved){
                LorVec_DiGamma = LorVec_gamma1 + LorVec_gamma2;
                LorVec_KaonPion = LorVec_K_plus + LorVec_DiGamma;
            }
        }

        //Recombination of Particles
        if(Kst2Kspiplus){
            LorVec_KshortMuMu = LorVec_K_short + LorVec_DiMuon;
        }
        if(Kst2Kpluspi0Merged){
            LorVec_KplusMuMu = LorVec_K_plus+ LorVec_DiMuon;
        }
        if(Kst2Kpluspi0Resolved){
            LorVec_KplusMuMu = LorVec_K_plus+ LorVec_DiMuon;
            LorVec_ReconBPiZero = LorVec_DiMuon + LorVec_K_plus + LorVec_pi_zero_resolved;
            dReconstructedBMFixPiZero = LorVec_ReconBPiZero.M();
            open_angle =LorVec_gamma1.Angle(LorVec_gamma2.Vect());
        }
        LorVec_K_star_plus = LorVec_KaonPion;
        LorVec_ReconB = LorVec_DiMuon + LorVec_K_star_plus;

        //double misIDed
        if(Kst2Kspiplus){
            LorVec_K_star_plus_From_DoubleMisIDed = LorVec_K_short + LorVec_mu_plus_MisIDed;
            LorVec_Jpsi_From_DoubleMisIDed = LorVec_pi_plus_MisIDed + LorVec_mu_minus;
            LorVec_B_plus_From_DoubleMisIDed = LorVec_Jpsi_From_DoubleMisIDed + LorVec_K_star_plus_From_DoubleMisIDed;
        }

        //save masses to doubles!
        dDiMuonM    = LorVec_DiMuon.M();
        dKaonPionM  = LorVec_KaonPion.M();
        dReconstructedBM = LorVec_ReconB.M();
        if(Kst2Kspiplus) dKshortMuMuM  = LorVec_KshortMuMu.M();


        //LL (Kshort decays inside the velo)
        if(Kst2Kspiplus){
            h_KshortDaughters->Fill(Ks_pi_minus_TRACK_Type, Ks_pi_plus_TRACK_Type);
            if(Ks_pi_minus_TRACK_Type == 3 && Ks_pi_plus_TRACK_Type == 3){
                nKshortDecayInVeLo = 1;
            }
            //DD (Kshort leaves velo)
            else nKshortDecayInVeLo = 0;

            if(nKshortDecayInVeLo)//cut on pion daughter PT for LL tracks!
              if(Ks_pi_plus_PT < 250. || Ks_pi_minus_PT < 250.)
                continue;
        }

        h_BplusMass_DiMuonMassRaw->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
        if(Kst2Kspiplus){
            if(nKshortDecayInVeLo == 0)h_BplusMass_DiMuonMassRawDD->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
            else h_BplusMass_DiMuonMassRawLL->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
        }

        //--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
        //trigger selection
        //L0 stage
        if(B_plus_L0MuonDecision_TOS == 0) continue;

        //HLT 1 stage
        if((B_plus_Hlt1TrackAllL0Decision_TOS == 0) && (B_plus_Hlt1TrackMuonDecision_TOS == 0)) continue;

        //HLT2 stage
        if((B_plus_Hlt2Topo2BodyBBDTDecision_TOS == 0) &&
                (B_plus_Hlt2Topo3BodyBBDTDecision_TOS == 0) &&
                (B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS == 0) &&
                (B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS == 0) &&
                //(B_plus_Hlt2SingleMuonDecision_TOS== 0) && //remove SingleMuon from selection, since its exclusive yield is 1% and pre-scaling messes up data/MC agreement
                (B_plus_Hlt2DiMuonDetachedDecision_TOS== 0) &&
                (B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS == 0) &&
                (B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS == 0)
                ) continue;

        nCounterTrigger++;

        //---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

        //some mass plots
        h_KstarMass2D->Fill(K_star_plus_M, LorVec_KaonPion.M());
        h_KstarMass->Fill(LorVec_KaonPion.M());

        if (Kst2Kpluspi0Merged){
            h_KGammaGammaMuMu_mass_before_preselection->Fill(dReconstructedBM);
            h_pi0_merged_candidates_pT->Fill(LorVec_pi_zero_merged.Perp());
            h_pi0_merged_candidates_M->Fill(LorVec_pi_zero_merged.M());

        }
        if (Kst2Kpluspi0Resolved){
            h_pi0_resolved_candidates_pT->Fill(LorVec_pi_zero_resolved.Perp());
            h_KGammaGammaMuMu_mass_before_preselection->Fill(dReconstructedBM);
            h_pi0_resolved_candidates_M->Fill(LorVec_DiGamma.M());
            h_pi0Mass2D->Fill(pi_zero_resolved_M,LorVec_DiGamma.M());
        }

        h_KstarMass2D->Fill(K_star_plus_M, LorVec_KaonPion.M());
        h_KstarMass->Fill(LorVec_KaonPion.M());

        //mass cuts
        if(TMath::Abs((UseDTF ? K_star_plus_M_DTF : K_star_plus_M) - PDGMASS.K_STAR_PLUS) > K_star_plus_mass_difference) continue; //remove non-resonant K+pi

        if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){

            if((UseDTF ? LorVec_K_star_plus.Pt() : K_star_plus_PT) < cut_K_star_plus_pT)continue;
            //if(TMath::Abs(LorVec_KaonPion.M() - PDGMASS.K_STAR_PLUS) > 150) continue; //remove non-resonant K+pi
        }
        h_KstarMasscut->Fill(LorVec_KaonPion.M());
        nCounterKstarMassWindow++;

        if((UseDTF ? B_plus_M_DTF : B_plus_M) < get_cut_B_plus_M_low(year) || (UseDTF ? B_plus_M_DTF : B_plus_M) > cut_B_plus_M_high)continue;
        //Check for non-DTF mass to allign with the older stripping
        if (Kst2Kpluspi0Resolved){
            if(B_plus_M < cut_B_plus_M_low_basic) continue;
        }

        //kinematic cuts
        if(LorVec_B_plus.Theta() <= cut_kin_Theta_low || LorVec_B_plus.Theta() >= cut_kin_Theta_up)continue;
        if(LorVec_K_star_plus.Theta() <= cut_kin_Theta_low || LorVec_K_star_plus.Theta() >= cut_kin_Theta_up)continue;
        if(LorVec_mu_minus.Theta() <= cut_kin_Theta_low || LorVec_mu_minus.Theta() >= cut_kin_Theta_up)continue;
        if(LorVec_mu_plus.Theta() <= cut_kin_Theta_low || LorVec_mu_plus.Theta() >= cut_kin_Theta_up)continue;

        /// newly added
        if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
            if ((UseDTF ? LorVec_B_plus.Pt() : B_plus_PT) < cut_B_plus_pT) continue;
            if (B_plus_DIRA_OWNPV < cut_B_plus_DIRA) continue;
        }

        if(Kst2Kspiplus){
            if(LorVec_K_short.Theta() <= cut_kin_Theta_low || LorVec_K_short.Theta() >= cut_kin_Theta_up)continue;
            if(LorVec_pi_plus.Theta() <= cut_kin_Theta_low || LorVec_pi_plus.Theta() >= cut_kin_Theta_up)continue;
            if(LorVec_KS_pi_plus.Theta() <= cut_kin_Theta_low || LorVec_KS_pi_plus.Theta() >= cut_kin_Theta_up)continue;
            if(LorVec_KS_pi_minus.Theta() <= cut_kin_Theta_low || LorVec_KS_pi_minus.Theta() >= cut_kin_Theta_up)continue;
        }
        if(Kst2Kpluspi0Merged){
            if(LorVec_K_plus.Theta() <= cut_kin_Theta_low || LorVec_K_plus.Theta() >= cut_kin_Theta_up)continue;
            if(LorVec_pi_zero_merged.Theta() <= cut_kin_Theta_low || LorVec_pi_zero_merged.Theta() >= cut_kin_Theta_up)continue;
            if((UseDTF ? LorVec_K_plus.Pt() : K_plus_PT) < 1350) continue;
        }
        if(Kst2Kpluspi0Resolved){
            if(LorVec_K_plus.Theta() <= cut_kin_Theta_low || LorVec_K_plus.Theta() >= cut_kin_Theta_up)continue;
            if(LorVec_pi_zero_resolved.Theta() <= cut_kin_Theta_low || LorVec_pi_zero_resolved.Theta() >= cut_kin_Theta_up)continue;
        }

        if(LorVec_mu_minus.Angle(LorVec_mu_plus.Vect()) <= cut_muon_angle)continue;

        if(Kst2Kspiplus){
//            if(LorVec_K_short.Angle(LorVec_pi_plus.Vect()) <= cut_muon_angle)continue;
            if(LorVec_pi_plus.Angle(LorVec_KS_pi_minus.Vect())    <= cut_muon_angle)continue;
            if(LorVec_pi_plus.Angle(LorVec_KS_pi_plus.Vect())     <= cut_muon_angle)continue;
            if(LorVec_KS_pi_plus.Angle(LorVec_KS_pi_minus.Vect()) <= cut_muon_angle)continue;
        }
        if(Kst2Kpluspi0Merged)if(LorVec_K_plus.Angle(LorVec_pi_zero_merged.Vect()) <= cut_muon_angle)continue;
        if(Kst2Kpluspi0Resolved)if(LorVec_K_plus.Angle(LorVec_pi_zero_resolved.Vect()) <= cut_muon_angle)continue;

        /// newly added
        if(Kst2Kpluspi0Merged)if((UseDTF ? LorVec_pi_zero_merged.Pt() : pi_zero_merged_PT) < cut_pi_zero_PT)continue;
        if(Kst2Kpluspi0Resolved)if((UseDTF ? LorVec_pi_zero_resolved.Pt() : pi_zero_resolved_PT) < cut_pi_zero_PT)continue;

        /// newly added
        //flight distance cuts
        if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
            if (B_plus_FDCHI2_OWNPV < cut_B_plus_FDCHI2) continue;
            if ((!Inc && !K1 && !B0 ) && K_star_plus_FDCHI2_OWNPV < cut_K_star_plus_FDCHI2) continue;
        }

        nCounterKinematic++;

        //primary vertex cut
        if (B_plus_IPCHI2_OWNPV      < cut_IPCHI2_OWNPV_low)continue;
        if (K_star_plus_IPCHI2_OWNPV < cut_IPCHI2_OWNPV_low)continue;
        if (mu_plus_IPCHI2_OWNPV     < cut_muon_IPCHI2_OWNPV_low)continue;
        if (mu_minus_IPCHI2_OWNPV    < cut_muon_IPCHI2_OWNPV_low)continue;

        if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged)if (B_plus_IPCHI2_OWNPV > cut_IPCHI2_OWNPV_high)continue;

        nCounterPrV++;

        /// newly added
       if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
            if (B_plus_ConePTasym < cut_B_plus_ConePTasym)continue;
       }

        //position vertex cuts?
        nCounterPoV++;

        h_BeforePIDcuts_ctl->Fill(TMath::Cos(ctl));
        h_BeforePIDcuts_ctk->Fill(TMath::Cos(ctk));
        h_BeforePIDcuts_phi->Fill(phi);

        //muon cuts
        if(!mu_plus_isMuon || !mu_minus_isMuon)continue;
        if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
            if (LorVec_DiMuon.M() > cut_DiMuon_M) continue; //Should be in stripping already
        }

        //particle ID
        if(mu_minus_ProbNNmu <= cut_muon_ProbNNmu)continue;
        if(mu_plus_ProbNNmu <= cut_muon_ProbNNmu)continue;
        if(Kst2Kspiplus){
            if(!pi_plus_hasRich)continue;
            //if(pi_plus_PIDK >= 25)continue;
        }

        if(Kst2Kpluspi0Merged){
            if(!K_plus_hasRich)continue;
            if(K_plus_ProbNNk <= cut_K_plus_ProbNNk)continue;
        }
        if(Kst2Kpluspi0Resolved){
            if(!K_plus_hasRich)continue;
            if(K_plus_ProbNNk <= cut_K_plus_ProbNNk)continue;
            if(gamma1_CL <= cut_gamma_CL)continue;
            if(gamma2_CL <= cut_gamma_CL)continue;
        }

        h_AfterPIDcuts_ctl->Fill(TMath::Cos(ctl));
        h_AfterPIDcuts_ctk->Fill(TMath::Cos(ctk));
        h_AfterPIDcuts_phi->Fill(phi);

        nCounterPID++;
        h_JpsiMass->Fill(LorVec_DiMuon.M());

        //only keep events inside q2 region
        if(LorVec_DiMuon.M2() < f_BinBoundaries[0] || LorVec_DiMuon.M2() > f_BinBoundaries[18])continue;

        nCounterDiMuonMassInQ2Range++;

        h_diMuonReconstructedBMass2D->Fill(dReconstructedBM,dDiMuonM);

        if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
            h_diMuonReconstructedBMass2DFixPiZero->Fill(dReconstructedBMFixPiZero,dDiMuonM);
        }

        if (Kst2Kpluspi0Merged){
            h_KGammaGammaMuMu_mass_after_preselection->Fill(dReconstructedBM);
        }
        if (Kst2Kpluspi0Resolved){
            h_KGammaGammaMuMu_mass_after_preselection->Fill(dReconstructedBM);
            h_KPiMuMu_KGammaGammaMuMu_mass->Fill(dReconstructedBMFixPiZero);
        }

        h_BplusMass_B02KsmumuMass->Fill(LorVec_KshortMuMu.M());

        if(Kst2Kspiplus){ //B0 veto
            if(TMath::Abs(LorVec_KshortMuMu.M() - PDGMASS.B_ZERO) < 50){
                //h_BplusMass_B02KsmumuMass->Fill(B_plus_M);
                nCounterB0KsmumuMass++;
                continue;
            }
        }

        ///////////////////////////////////////
        ///////////////////////////////////////
        ///////////////////////////////////////

        nDiMuonMassBin = 0;
        Is_Phi_JPsi_Psi2S = 0;
        while(LorVec_DiMuon.M2() >= f_BinBoundaries[nDiMuonMassBin]){
            nDiMuonMassBin++;
        }

        if(nDiMuonMassBin == 2){
            Is_Phi_JPsi_Psi2S = 1;
        }
        else if(nDiMuonMassBin == 9){
            h_JpsiMassResonance->Fill(LorVec_DiMuon.M());
            h_CharmoniumResonance->Fill(LorVec_DiMuon.M());
            h_ReconstructBmassInCharmoniumResonance->Fill(LorVec_ReconB.M());
            Is_Phi_JPsi_Psi2S = 1;
        }
        else if(nDiMuonMassBin == 12){
            h_Psi2sMassResonance->Fill(LorVec_DiMuon.M());
            h_CharmoniumResonance->Fill(LorVec_DiMuon.M());
            h_ReconstructBmassInCharmoniumResonance->Fill(LorVec_ReconB.M());
            Is_Phi_JPsi_Psi2S = 1;
        }
        else{
            h_JpsiMassCut->Fill(LorVec_DiMuon.M());
            h_ReconstructBmassOffResonance->Fill(LorVec_ReconB.M());
            nCounterJpsiResonance++;
        }

        /*if((LorVec_DiMuon.M2() < 8e6 || LorVec_DiMuon.M2() > 11e6) && (LorVec_DiMuon.M2() < 12.5e6 || LorVec_DiMuon.M2() > 15e6) && (LorVec_DiMuon.M2() < 0.98e6 || LorVec_DiMuon.M2() > 1.1e6)){
            h_JpsiMassCut->Fill(LorVec_DiMuon.M());
            nCounterJpsiResonance++;
        }*/

        ///////////////////////////////////////
        if(Kst2Kpluspi0Merged || Kst2Kpluspi0Resolved){
            // if (nDiMuonMassBin == 9 || nDiMuonMassBin == 12) continue;
            nCounterB0KsmumuMass++;
            //h_BplusMass_B02KsmumuMass->Fill(B_plus_M);
        }

        ///very soft cuts to ensure a more even distribution for the BDT training:
        if(Kst2Kspiplus){
          if((UseDTF ? LorVec_B_plus.Pt() : B_plus_PT) > 55000.)continue;
          if((UseDTF ? LorVec_K_short.Pt() : K_short_PT) > (nKshortDecayInVeLo == 1 ? 9000. : 14000.))continue;
          if(TMath::Log(B_plus_IPCHI2_OWNPV) < -8.)continue;
        }

        //get subset number from modulus(3) of the eventnumber
        //range 0 to 2
        nSubSetNumber = nEventNumber % 3;

        //plot the double misIDed reconstructed mass distributions for Kshort Channel:
        if(Kst2Kspiplus){

            dJpsiMisIDed_M = LorVec_Jpsi_From_DoubleMisIDed.M();

            h_Jpsi_Spectrum_MisIDed_Muon->Fill(dJpsiMisIDed_M);
            h_BplusM_vs_BplusM_DoubleMisIDed->Fill(LorVec_B_plus_From_DoubleMisIDed.M(), dReconstructedBM);
            h_BplusM_vs_JpsiM_DoubleMisIDed->Fill(LorVec_B_plus_From_DoubleMisIDed.M(), dJpsiMisIDed_M);
            h_KstarM_vs_JpsiM_DoubleMisIDed->Fill(LorVec_K_star_plus_From_DoubleMisIDed.M(), dJpsiMisIDed_M);
            if(nDiMuonMassBin != 2 && nDiMuonMassBin != 9 && nDiMuonMassBin != 12){
                h_Q2_Spectrum_MisIDed_Muon->Fill(LorVec_Jpsi_From_DoubleMisIDed.M2() * 1e-6);
                h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->Fill(dJpsiMisIDed_M);
            }
        }

        //remove double mis-identified events:
        if(Kst2Kspiplus){
            if(TMath::Abs(dJpsiMisIDed_M - PDGMASS.J_PSI) < 50 && TMath::Abs(LorVec_K_star_plus_From_DoubleMisIDed.M() - PDGMASS.K_STAR_PLUS) < 30){
                n_DoubleMisID++;
                continue;
            }
        }

        //plot for 2D histos
        h_BplusMass_DiMuonMass->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
        if(Kst2Kspiplus){
            if(nKshortDecayInVeLo == 1)h_BplusMass_DiMuonMassLL->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
            else h_BplusMass_DiMuonMassDD->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
        }

        dQ2 = LorVec_DiMuon.M2();
        //---------------------------------------------------------------------------------

        summary_tree->Fill();
    }

    output->cd();

    coutInfo("Remove " + to_string(n_DoubleMisID)+ " double mis-IDed events");
    coutInfo("Save " + to_string(summary_tree->GetEntries())+ " events to new file!");

    summary_tree->Write("",TObject::kWriteDelete);
    output->Close();

    coutInfo("Save quality histograms to seperate file!");

    histos->cd();
    h_JpsiMass->Write("",TObject::kWriteDelete);
    h_JpsiMassCut->Write("",TObject::kWriteDelete);
    h_JpsiMassResonance->Write("",TObject::kWriteDelete);
    h_Psi2sMassResonance->Write("",TObject::kWriteDelete);
    h_CharmoniumResonance->Write("",TObject::kWriteDelete);

    h_ReconstructBmassInCharmoniumResonance->Write("",TObject::kWriteDelete);
    h_ReconstructBmassOffResonance->Write("",TObject::kWriteDelete);

    h_KstarMass->Write("",TObject::kWriteDelete);
    h_KstarMasscut->Write("",TObject::kWriteDelete);
    h_KstarMass2D->Write("",TObject::kWriteDelete);

    h_BplusMass_B02KsmumuMass->Write("",TObject::kWriteDelete);

    h_KshortDaughters->Write("",TObject::kWriteDelete);
    h_Jpsi_Spectrum_MisIDed_Muon->Write("",TObject::kWriteDelete);
    h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->Write("",TObject::kWriteDelete);
    h_Q2_Spectrum_MisIDed_Muon->Write("",TObject::kWriteDelete);
    h_BplusM_vs_BplusM_DoubleMisIDed->Write("",TObject::kWriteDelete);
    h_BplusM_vs_JpsiM_DoubleMisIDed->Write("",TObject::kWriteDelete);
    h_KstarM_vs_JpsiM_DoubleMisIDed->Write("",TObject::kWriteDelete);

    h_pi0_merged_candidates_M->Write("",TObject::kWriteDelete);
    h_pi0_resolved_candidates_M->Write("",TObject::kWriteDelete);
    h_pi0_merged_candidates_pT->Write("",TObject::kWriteDelete);
    h_pi0_resolved_candidates_pT->Write("",TObject::kWriteDelete);
    h_KGammaGammaMuMu_mass_before_preselection->Write("",TObject::kWriteDelete);
    h_KGammaGammaMuMu_mass_after_preselection->Write("",TObject::kWriteDelete);
    h_KPiMuMu_KGammaGammaMuMu_mass->Write("",TObject::kWriteDelete);
    h_pi0Mass2D->Write("",TObject::kWriteDelete);
    h_diMuonReconstructedBMass2D->Write("",TObject::kWriteDelete);
    h_diMuonReconstructedBMass2DFixPiZero->Write("",TObject::kWriteDelete);
    
    h_BplusMass_DiMuonMass->Write("",TObject::kWriteDelete);
    h_BplusMass_DiMuonMassDD->Write("",TObject::kWriteDelete);
    h_BplusMass_DiMuonMassLL->Write("",TObject::kWriteDelete);
    h_BplusMass_DiMuonMassRaw->Write("",TObject::kWriteDelete);
    h_BplusMass_DiMuonMassRawDD->Write("",TObject::kWriteDelete);
    h_BplusMass_DiMuonMassRawLL->Write("",TObject::kWriteDelete);

    h_BeforePIDcuts_ctl->Write("",TObject::kWriteDelete);
    h_BeforePIDcuts_ctk->Write("",TObject::kWriteDelete);
    h_BeforePIDcuts_phi->Write("",TObject::kWriteDelete);
    h_AfterPIDcuts_ctl->Write("",TObject::kWriteDelete);
    h_AfterPIDcuts_ctk->Write("",TObject::kWriteDelete);
    h_AfterPIDcuts_phi->Write("",TObject::kWriteDelete);

    histos->Close();

    coutInfo("All done for " + year + " " + magnet + string(MC ? (PHSP ? " (PHSP)" : (ReferenceChannel ? " (RefMC) " : " (MC)") ) : "") + "!");

    delete output;
    delete histos;

    //Add TM branches for TMVA running
    if(!MC) copyDataForTMVA(year,magnet);

    return 1;

}

//run multiple or all functions:
int runAllYearsAndMagPol(bool MC = false, bool ReferenceChannel = false, bool PHSP = false,  bool B0 = false,  bool K1 = false, bool Inc = false, Int_t Run = 1){

    std::vector<string> years = yearsVector(MC,ReferenceChannel,PHSP,Run);

    string Magnets[] = {"down", "up"};

    for(unsigned int y = 0; y < years.size(); y++){ //loop over all years
      for(string & magnet : Magnets){ //loop over both magnet polaritys
        if(preSelection(years.at(y), magnet, MC, ReferenceChannel, PHSP, B0, K1, Inc) == 0){
          coutERROR("Failed for " + years.at(y) + " " + magnet);
          return 0;
        }
      }
    }
    return 1;

}

//signal channels
int runAllSignalData(Int_t Run = 1){
    if(runAllYearsAndMagPol(false, false, false, false, false, false, Run) == 0){
        coutERROR("Failed for signal data sample of decay " + TheDecay + "!");
        return 0;
    }
    coutInfo("All signal data for decay " + TheDecay + " successfully processed!!");
    return 1;
}

int runAllSignalMC(Int_t Run = 1){
    if(runAllYearsAndMagPol(true, false, false, false, false, false, Run) == 0){
        coutERROR("Failed for signal MC sample of decay " + TheDecay + "!");
        return 0;
    }
    coutInfo("All signal MC for decay " + TheDecay + " successfully processed!!");
    return 1;
}

//reference MC channels
int runAllRefMC(Int_t Run = 1){
    if(runAllYearsAndMagPol(true, true, false, false, false, false, Run) == 0){
        coutERROR("Failed for reference MC sample of decay " + TheDecay + "!");
        return 0;
    }
    coutInfo("All reference MC for decay " + TheDecay + " successfully processed!!");
    return 1;
}

//reference MC channels
int runAllPHSPMC(Int_t Run = 1){
    if(runAllYearsAndMagPol(true, false, true, false, false, false, Run) == 0){
        coutERROR("Failed for PSHP MC sample of decay " + TheDecay + "!" );
        return 0;
    }
    coutInfo("All PHSP MC for decay " + TheDecay +" successfully processed!!" );
    return 1;
}

//B0toKst MC channels
int runAllB0MC(Int_t Run = 1, bool ReferenceChannel = false){
        //BtoXMuMu not available at the moment, set ReferenceChannel to true
    if(runAllYearsAndMagPol(true, ReferenceChannel, false, true, false, false, Run) == 0){
        coutERROR("Failed for inclusive MC sample of decay " + TheDecay + "!" );
       // return 0;
    }
    coutInfo("All inclusive MC for decay " + TheDecay +" successfully processed!!" );
    return 1;
}
//BtoK1mumu MC channels
int runAllK1MC(bool ReferenceChannel = false){
   if(preSelection("2012" , "down", true, ReferenceChannel, false, false, true, false) == 0){
       coutERROR("Failed for K1 MC sample of decay " + TheDecay + "!" );
        return 0;  // Easier than specifying only 2012 available
    }
    if(preSelection("2012" , "up", true, ReferenceChannel, false, false, true, false) == 0){
       coutERROR("Failed for K1 MC sample of decay " + TheDecay + "!" );
        return 0;  // Easier than specifying only 2012 available
    }
    coutInfo("All K1 MC for decay " + TheDecay +" successfully processed!!" );
    return 1;
}
//Inclusive MC channels
int runAllIncMC(Int_t Run = 1){

    //BtoXMuMu not available at the moment, set ReferenceChannel to true
    bool ReferenceChannel = true;

    string Magnets[] = {"down", "up"};

    for(auto yr: yearsInc(Run)){ //loop over all years
      for(string & magnet : Magnets){ //loop over both magnet polaritys
        if(preSelection(yr, magnet, true, ReferenceChannel, false, false, false, true) == 0){
          coutERROR("Failed for " + yr + " " + magnet);
          return 0;
        }
      }
    }
    coutInfo("All inclusive MC for decay " + TheDecay +" successfully processed!!" );
    return 1;
}

//All background MC channels
int runAllBkgMC(){
    runAllB0MC(1,0);
    runAllB0MC(1,1);
    runAllB0MC(2,0);
    runAllB0MC(2,1);

    runAllK1MC(0); //only 2012 available
    runAllK1MC(1); //only 2012 available

    runAllIncMC(12);

    coutInfo("All inclusive MC for decay " + TheDecay +" successfully processed!!" );
    return 1;
}


int getTCKs(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false) {

    gStyle -> SetOptStat(0);
    LHCbStyle();
    gROOT->SetBatch(kTRUE);
    
    TChain* tree = nullptr;

    if(ReferenceChannel)MC = true;
    if(PHSP)MC = true;

    if(ReferenceChannel && PHSP){
        coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
        PHSP = false;
    }

    //Kst2Kpluspi0Resolved case
    if(Kst2Kpluspi0Resolved){
        tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    //Kst2Kpluspi0Merged case
    if(Kst2Kpluspi0Merged){
        tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    //Kst2Kspiplus case
    if(Kst2Kspiplus){
        tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    if(tree == nullptr){
        coutERROR("No files found!");
        return 0;
    }


    int N = tree->GetEntries();
    if(N == 0)coutERROR("No events found in file!");
    else coutInfo("Old files contain " + to_string(N)+ " events");

    //Disable all branches
    tree->SetBranchStatus("*",0);

    //activate needed branch
    UInt_t HLT1TCK;
    tree->SetBranchStatus("HLT1TCK",1);
    tree->SetBranchAddress("HLT1TCK", &HLT1TCK);

    coutInfo("Get Min and Max from HLT1TCK branch");
    Double_t HLT1TCKMin = tree->GetMinimum("HLT1TCK");
    Double_t HLT1TCKMax = tree->GetMaximum("HLT1TCK");

    UInt_t runNumber;
    tree->SetBranchStatus("runNumber",1);
    tree->SetBranchAddress("runNumber", &runNumber);

    coutInfo("Get Min and Max from runNumber branch");
    Double_t runNumberMin = tree->GetMinimum("runNumber");
    Double_t runNumberMax = tree->GetMaximum("runNumber");

    TCanvas * cRunNumberTCKs = new TCanvas("RunNumberTCKs", "RunNumberTCKs");
    cRunNumberTCKs->cd();
    TH2D * hRunNumberTCK = new TH2D(Form("%s_RunNumber_vs_TCKs_used_%s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()),Form("%s: RunNumber vs. TCKs used in %s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()), 100, runNumberMin, runNumberMax, 100, HLT1TCKMin, HLT1TCKMax);

    //vector to count total number of each TCK. Index 0 is the sum of all TCKs, so should be equal to number of events
    std::vector<Int_t>TCKcounts;
    std::vector<UInt_t>FoundTCKs;

    coutInfo("Start loop over " + to_string(N) + " events!");

    //Fill in the first entry
    tree->GetEntry(0);
    hRunNumberTCK->Fill(runNumber, HLT1TCK);
    FoundTCKs.push_back(HLT1TCK);
    TCKcounts.push_back(1);

    //Check for different TCKs in a loop
    for(int ev = 1; ev < N; ev++){

        tree->GetEntry(ev);
        hRunNumberTCK->Fill(runNumber, HLT1TCK);

        if(ev%100000 == 0)std::cout << Form("[COUNTER]\t\t %i / %i", ev, N) << std::endl;

        unsigned int TCKindex = 0;

        //If the TCK is different to all previous TCKs, fill it into FoundTCKs
        while(FoundTCKs.at(TCKindex) != HLT1TCK){
            TCKindex++;
            if(FoundTCKs.size() == TCKindex){
                FoundTCKs.push_back(HLT1TCK);
                TCKcounts.push_back(0);
                TCKindex--;
                break;
            }
        }
        //Add count to a given TCK
        TCKcounts.at(TCKindex)++;

    }

    if(TCKcounts.size() != FoundTCKs.size()){
        coutERROR("Number of found TCKs and number of counters is not equal! Exit program...");
        return 0;
    }

    unsigned int i = 1;//, j = 0;
    Int_t TmpTCKcount, TmpFoundTCK;

    // order TCKs in order (says David)
    while(i < TCKcounts.size()){

        unsigned int j = i;

        while(FoundTCKs.at(j) < FoundTCKs.at(j - 1)){

            TmpFoundTCK = FoundTCKs.at(j);
            FoundTCKs.at(j) = FoundTCKs.at(j - 1);
            FoundTCKs.at(j - 1) = TmpFoundTCK;

            TmpTCKcount = TCKcounts.at(j);
            TCKcounts.at(j) = TCKcounts.at(j - 1);
            TCKcounts.at(j - 1) = TmpTCKcount;

            j--;

            if(j ==0)break;
        }

        i++;
    }

    bool print_latex = true;
    if(print_latex)
        std::cout << year << " " << magnet;
    else{
        coutInfo("Number of events of each used TCKs in " + year + " " + magnet +":");
        std::cout << "[INFO]" << std::endl;
        coutInfo("TCK\tEvtCounter\tPercentage");
        std::cout << "[INFO]" << std::endl;
    }

    for(unsigned int index = 0; index < TCKcounts.size(); index++){
        //std::cout << Form("[INFO]\t\t%d \t\t%d\t\t%.2f", FoundTCKs.at(index), TCKcounts.at(index), 100./N*TCKcounts.at(index)) << std::endl;
        if(print_latex)
            std::cout << "\t&" << hex << FoundTCKs.at(index) << "\t&" << dec << TCKcounts.at(index) << "\t&" << std::setprecision(2) << std::fixed << 100./N*TCKcounts.at(index) << "\\\\" << std::endl;
        else
            std::cout << "[INFO]" << hex << FoundTCKs.at(index) << "\t" << dec << TCKcounts.at(index) << "\t" << 100./N*TCKcounts.at(index) << std::endl;
    }
    std::cout << "[INFO]" << std::endl;

    TCanvas * cTCKs = new TCanvas("TCKs", "TCKs");
    TH1D * hTCKs = new TH1D(Form("%s_TCKs_used_%s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()),Form("%s: TCKs used in %s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()), TCKcounts.size(), 0, TCKcounts.size());

    //Label TCK bins with TCK
    for(unsigned bin = 0; bin < TCKcounts.size(); bin++){
        for(int counts = 0; counts < TCKcounts.at(bin); counts++){
            hTCKs->Fill(bin);
        }

        hTCKs->GetXaxis()->SetBinLabel(bin + 1, Form("%x", FoundTCKs.at(bin)));
        //hTCKs->GetXaxis()->ChangeLabel(bin + 1, 45., -1., -1, -1, -1, Form("%i", FoundTCKs.at(bin)));

    }

    //Save the TCK histograms
    gStyle->SetOptStat(0);

    string suffix = "";
    if(MC){
        if(PHSP)suffix.append("_PHSP");
        else if(ReferenceChannel)suffix.append("_Ref");
        else suffix.append("_MC");
    }

    cTCKs->cd();
    hTCKs->GetXaxis()->SetTitle("TCKid");
    hTCKs->GetYaxis()->SetTitle("Events / TCK");
    hTCKs->SetMarkerStyle(kFullCircle);
    hTCKs->Draw("PE");
    cTCKs->SaveAs(Form("%s/TCKinfo/TCKs%s%s%s%s.eps",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));
    cTCKs->SaveAs(Form("%s/TCKinfo/TCKs%s%s%s%s.root",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));

    cRunNumberTCKs->cd();
    hRunNumberTCK->GetYaxis()->SetTitle("TCKid");
    hRunNumberTCK->GetXaxis()->SetTitle("RunNumber");
    hRunNumberTCK->Draw("COLZ");
    cRunNumberTCKs->SaveAs(Form("%s/TCKinfo/RunNumberTCKs%s%s%s%s.eps",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));
    cRunNumberTCKs->SaveAs(Form("%s/TCKinfo/RunNumberTCKs%s%s%s%s.root",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));

    return 1;

}

int ExclusiveTriggerLineSelection(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false) {

    gStyle -> SetOptStat(0);
    LHCbStyle();
    gROOT->SetBatch(kTRUE);

    TChain* tree = nullptr;

    if(ReferenceChannel)MC = true;
    if(PHSP)MC = true;

    if(ReferenceChannel && PHSP){
        coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
        PHSP = false;
    }

    //Kst2Kpluspi0Resolved case
    if(Kst2Kpluspi0Resolved){
        tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    //Kst2Kpluspi0Merged case
    if(Kst2Kpluspi0Merged){
        tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    //Kst2Kspiplus case
    if(Kst2Kspiplus){
        tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
        tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
    }

    if(tree == nullptr){
        coutERROR("No files found!");
        return 0;
    }


    int N = tree->GetEntries();
    if(N == 0)coutERROR("No events found in file!");
    else coutInfo("Old files contain " + to_string(N)+ " events");

    //Disable all branches
    tree->SetBranchStatus("*",0);

    //activate needed branch
    //L0
    Bool_t B_plus_L0MuonDecision_TOS = false;

    //HLT1
    Bool_t B_plus_Hlt1TrackAllL0Decision_TOS = false;
    Bool_t B_plus_Hlt1TrackMuonDecision_TOS = false;

    //HLT2
    Bool_t B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS = false;
    Bool_t B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS = false;

    Bool_t B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS = false;
    Bool_t B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS = false;

    Bool_t B_plus_Hlt2Topo2BodyBBDTDecision_TOS = false;
    Bool_t B_plus_Hlt2Topo3BodyBBDTDecision_TOS = false;

    Bool_t B_plus_Hlt2SingleMuonDecision_TOS = false;

    Bool_t B_plus_Hlt2DiMuonDetachedDecision_TOS = false;


    tree -> SetBranchAddress( "B_plus_L0MuonDecision_TOS" , &B_plus_L0MuonDecision_TOS );

    if(year == "2011" || year == "2012"){
        tree -> SetBranchAddress( "B_plus_Hlt1TrackAllL0Decision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );
    }
    else if(year == "2015" || year == "2016" || year == "2017" || year == "2018"){
        tree -> SetBranchAddress( "B_plus_Hlt1TrackMVADecision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );

        tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
        tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );

        if(year == "2016" || year == "2017" || year == "2018"){ //found to be FALSE for all 2015 data
            tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS );
            tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS );
        }
    }
    else{
        coutERROR("Incorrect year given for the assosiation of HLT1 lines with branches! Year given: " + year);
        return 0;
    }
    tree -> SetBranchAddress( "B_plus_Hlt1TrackMuonDecision_TOS" , &B_plus_Hlt1TrackMuonDecision_TOS );
    tree -> SetBranchAddress( "B_plus_Hlt2SingleMuonDecision_TOS" , &B_plus_Hlt2SingleMuonDecision_TOS );
    tree -> SetBranchAddress( "B_plus_Hlt2DiMuonDetachedDecision_TOS" , &B_plus_Hlt2DiMuonDetachedDecision_TOS );

    //count exclusive trigger decisions:

    Int_t ex_Hlt1TrackAllL0 = 0,
          ex_Hlt1TrackMuon = 0;
    Int_t ex_Hlt2SingleMuon = 0,
          ex_Hlt2DiMuonDetached = 0,
          ex_Hlt2TopoMu2Body = 0,
          ex_Hlt2TopoMu3Body = 0,
          ex_Hlt2Topo2Body = 0,
          ex_Hlt2Topo3Body = 0,
          ex_Hlt2TopoMuMu2Body = 0,
          ex_Hlt2TopoMuMu3Body = 0;

    coutInfo("Start loop over " +to_string(N)+ " events!");

    for(int ev = 0; ev < N; ev++){
        tree->GetEntry(ev);
        if(ev%100000 == 0)std::cout << Form("[COUNTER]\t\t %i / %i", ev, N) << std::endl;

        //HLT1
        if(B_plus_Hlt1TrackAllL0Decision_TOS && !B_plus_Hlt1TrackMuonDecision_TOS)
            ex_Hlt1TrackAllL0++;
        if(!B_plus_Hlt1TrackAllL0Decision_TOS && B_plus_Hlt1TrackMuonDecision_TOS)
            ex_Hlt1TrackMuon++;

        //HLT2
        if(B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
            ex_Hlt2SingleMuon++;
        else if(!B_plus_Hlt2SingleMuonDecision_TOS && B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
            ex_Hlt2DiMuonDetached++;
        else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
            ex_Hlt2TopoMu2Body++;
        else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
            ex_Hlt2TopoMu3Body++;
        else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
            ex_Hlt2Topo2Body++;
        else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
            ex_Hlt2Topo3Body++;
        else if(year == "2016" || year == "2017" || year == "2018"){
            if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
                ex_Hlt2TopoMuMu2Body++;
            else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
                ex_Hlt2TopoMuMu3Body++;
        }
    }

    bool print_latex = true;

    std::cout << std::endl;
    std::cout << year << " " << magnet << ":" << std::endl;
    std::cout << "Trigger line\texclusive events\tpercentage" << std::endl;
    std::cout << "HLT1" << std::endl;
    std::cout << "Hlt1TrackAllL0\t\t" << (print_latex ? " &" : "") << ex_Hlt1TrackAllL0 << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt1TrackAllL0 / N << std::endl;
    std::cout << "Hlt1TrackMuon\t\t" << (print_latex ? " &" : "") << ex_Hlt1TrackMuon << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt1TrackMuon / N << std::endl;
    std::cout << "HLT2" << std::endl;
    std::cout << "Hlt2SingleMuon\t\t" << (print_latex ? " &" : "") << ex_Hlt2SingleMuon << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2SingleMuon / N << std::endl;
    std::cout << "Hlt2DiMuonDetached\t" << (print_latex ? " &" : "") << ex_Hlt2DiMuonDetached << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2DiMuonDetached / N << std::endl;
    std::cout << "Hlt2TopoMu2Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMu2Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMu2Body / N << std::endl;
    std::cout << "Hlt2TopoMu3Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMu3Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMu3Body / N << std::endl;
    std::cout << "Hlt2Topo2Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2Topo2Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2Topo2Body / N << std::endl;
    std::cout << "Hlt2Topo3Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2Topo3Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2Topo3Body / N << std::endl;
    std::cout << "Hlt2TopoMuMu2Body\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMuMu2Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMuMu2Body / N << std::endl;
    std::cout << "Hlt2TopoMuMu3Body\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMuMu3Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMuMu3Body / N << std::endl;

    return 1;

}

Int_t plotHistos(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false){

    //Reference channel and PHSP are both MC data
    if(ReferenceChannel)MC = true;
    if(PHSP)MC = true;
    
    if(ReferenceChannel && PHSP){
        coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
        PHSP = false;
    }
    

    gStyle -> SetOptStat(0);
    LHCbStyle();

    TFile * histos = nullptr;
    if(Kst2Kpluspi0Resolved){
        if(!MC)histos = new TFile(Form("%s/data/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
        else{
        if(ReferenceChannel)histos = new TFile(Form("%s/data/MC/RefKplusPi0Resolved/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
            else if(PHSP)histos = new TFile(Form("%s/data/PHSP/KplusPi0Resolved/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
            else histos = new TFile(Form("%s/data/MC/KplusPi0Resolved/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
        }
    }
    if(Kst2Kspiplus){
        if(!MC)histos = new TFile(Form("%s/data/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
        else{
            if(ReferenceChannel)histos = new TFile(Form("%s/data/MC/RefKshortPiplus/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
	    else if(PHSP)histos = new TFile(Form("%s/data/PHSP/KshortPiplus/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
            else histos = new TFile(Form("%s/data/MC/KshortPiplus/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
        }
    }

    if(histos == nullptr){
    coutERROR("File with histograms not found or opened! Exit...");
	return 0;
    }

    TCanvas * c1 = new TCanvas("c1", "c1");
    c1->cd();
    TH2D * h_KstarM_vs_JpsiM_DoubleMisIDed = nullptr;
    h_KstarM_vs_JpsiM_DoubleMisIDed = (TH2D * ) histos->Get("h_KstarM_vs_JpsiM_DoubleMisIDed");
    if(h_KstarM_vs_JpsiM_DoubleMisIDed == nullptr){
    coutERROR("Histogram 'h_KstarM_vs_JpsiM_DoubleMisIDed' not found file! Exit...");
        return 0;
    }
    h_KstarM_vs_JpsiM_DoubleMisIDed->GetXaxis()->SetTitle((UseDTF ? "m([#mu^{+}]_{#pi^{+}}K_{S}^{0})" : "m([#mu^{+}]_{#pi^{+}}#pi^{-}#pi^{+})"));
    h_KstarM_vs_JpsiM_DoubleMisIDed->GetYaxis()->SetTitle("m([#pi^{+}]_{#mu^{+}}#mu^{-})");

    h_KstarM_vs_JpsiM_DoubleMisIDed->Draw("COLZ");

    //insert veto-region as overlay:
    Float_t xLow  = PDGMASS.K_STAR_PLUS - 30.;
    Float_t xHigh = PDGMASS.K_STAR_PLUS + 30.;
    Float_t yLow  = PDGMASS.J_PSI - 50.;
    Float_t yHigh = PDGMASS.J_PSI + 50.;

    TGraphErrors * box = new TGraphErrors();
    
    box->SetPoint(0, PDGMASS.K_STAR_PLUS, PDGMASS.J_PSI);
    box->SetPointError(0, 30., 50.); 
    box->SetFillColor(kRed);
    box->SetFillStyle(3445);
    box->Draw("2SAME"); 

    TLine * boxframe[4];
    boxframe[0] = new TLine(xLow,  yLow,  xHigh, yLow);
    boxframe[1] = new TLine(xHigh, yLow,  xHigh, yHigh);
    boxframe[2] = new TLine(xHigh, yHigh, xLow,  yHigh);
    boxframe[3] = new TLine(xLow,  yHigh, xLow,  yLow); 

    for(int i = 0; i < 4; i++){
	boxframe[i]->SetLineWidth(2);
	boxframe[i]->SetLineColor(kRed);
	boxframe[i]->Draw("SAME");
    }

    c1->Print(Form("%s/PreSelectionPlots/%s_%s%s%s%s%s_DoubleMisID_2Dplot.eps",
    thePath.c_str(),
	TheDecay.c_str(),
	year.c_str(),
	magnet.c_str(), 
	(MC ? "_MC" : ""), 
	(ReferenceChannel ? "_Ref" : ""),
	(PHSP ? "_PHSP" : "")));

    TCanvas * c2 = new TCanvas("c2", "c2");
    c2->cd();
    TH1D * h_BplusMass_B02KsmumuMass = nullptr;
    h_BplusMass_B02KsmumuMass = (TH1D * ) histos->Get("h_BplusMass_B02KsmumuMass");
    if(h_BplusMass_B02KsmumuMass == nullptr){
        coutERROR("Histogram 'h_BplusMass_B02KsmumuMass' not found file! Exit...");
        return 0;
    }
    h_BplusMass_B02KsmumuMass->GetXaxis()->SetTitle((UseDTF ? "m(K_{S}^{0}#mu^{+}#mu^{-})" : "m(#pi^{-}#pi^{+}#mu^{+}#mu^{-})"));
    h_BplusMass_B02KsmumuMass->GetYaxis()->SetTitle("candidates / 20MeV");

    Float_t MaxContent = -1.;
    for(int b = 1; b <= h_BplusMass_B02KsmumuMass->GetXaxis()->GetNbins(); b++){
	if(h_BplusMass_B02KsmumuMass->GetBinContent(b) > MaxContent)MaxContent = h_BplusMass_B02KsmumuMass->GetBinContent(b);
    }
    h_BplusMass_B02KsmumuMass->GetYaxis()->SetRangeUser(0, MaxContent * 1.2);

    h_BplusMass_B02KsmumuMass->Draw();

    TGraphErrors * box2 = new TGraphErrors();

    box2->SetPoint(0, PDGMASS.B_ZERO, MaxContent * 0.6);
    box2->SetPointError(0, 50., MaxContent * 0.6);
    box2->SetFillColor(kRed);
    box2->SetFillStyle(3445);
    box2->Draw("2SAME");

    TLine * vetoframe[2];
    vetoframe[0] = new TLine(PDGMASS.B_ZERO - 50., 0., PDGMASS.B_ZERO - 50., MaxContent * 1.2);
    vetoframe[1] = new TLine(PDGMASS.B_ZERO + 50., 0., PDGMASS.B_ZERO + 50., MaxContent * 1.2);
    for(int i = 0; i < 2; i++){
        vetoframe[i]->SetLineWidth(2);
        vetoframe[i]->SetLineColor(kRed);
        vetoframe[i]->Draw("SAME");
    }


    c2->Print(Form("%s/PreSelectionPlots/%s_%s%s%s%s%s_B0plusRandomPion.eps",
        thePath.c_str(),
        TheDecay.c_str(),
        year.c_str(),
        magnet.c_str(),
        (MC ? "_MC" : ""),
        (ReferenceChannel ? "_Ref" : ""),
        (PHSP ? "_PHSP" : "")));


    return 1;

}

void compare_PIDgen_Eventnumbers(int year = 2011, string magnet = "down"){

    string subDecay = Kst2Kspiplus ? "KshortPiplus" : "KplusPi0Resolved";

    TChain * treePIDgen, * tree;
    if(Kst2Kspiplus){
        treePIDgen = new TChain("b2KstKs0pimumu_Tuple/DecayTree");
        tree       = new TChain("b2KstKs0pimumu_Tuple/DecayTree");
    }
    else{
        treePIDgen = new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
        tree       = new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
    }

    treePIDgen->Add((path_to_data + "/data/PHSP/" + subDecay + "/" + std::to_string(year) + magnet + "/*B2KstmumuPIDgen_*.root").c_str());
    tree      ->Add((path_to_data + "/data/PHSP/" + subDecay + "/" + std::to_string(year) + magnet + "/*B2Kstmumu_*.root").c_str());

    std::cout << "------------------------------------" << std::endl;
    std::cout << "----" << year << " " << magnet << "----" << std::endl;
    std::cout << "------------------------------------" << std::endl;
    std::cout << "PIDgen tuple entries=" << treePIDgen->GetEntries() << std::endl;
    std::cout << "Normal tuple entries=" << tree      ->GetEntries() << std::endl;
    std::cout << "------------------------------------" << std::endl; 
}
