#include "../GlobalFunctions.hh"
#include "../Paths.hpp"
#include <algorithm> //Finding stuff in a vector

//////////////////////////////////////////////////////////////////////
/// sWeightSelection()
/// loops over all .root-files created by BDT selection. More kinematic, PID and
/// some more cuts are applied and saved in new file.
///

int sWeightSelection(bool useOld = false, std::string year = "2011", std::string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false) {

    TFile* input = 0;
    TFile* output_new = 0;
    TFile* output_old = 0;
    TTree *tree;

    if (PHSP) MC = true; //phase-space is always MC

    string input_path = GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,false);
    coutDebug("input_path: " + input_path);

    string output_new_path = GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,false);
    replace(output_new_path,".root","_tmp.root");
    string output_old_path = GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,false);
    replace(output_old_path,".root","_old.root");



    input       = new TFile(input_path.c_str(),"UPDATE");
    output_new  = new TFile(output_new_path.c_str(),"RECREATE");
    output_old  = new TFile(output_old_path.c_str(),useOld ? "UPDATE" : "RECREATE");


    if (useOld) tree = (TTree*)output_old->Get("DecayTree");
    else        tree = (TTree*)input->Get("DecayTree");

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

        //copy old tree into a backup
    output_old->cd();
    TTree *cpTree = tree->CloneTree();
    if (!cpTree) coutERROR("Tree not cloned!");
    else coutDebug("Tree copied!");
    cpTree->Write("",TObject::kOverwrite); //copy old trees into a backup; prevent duplicate writting
    if (!useOld) output_old->Close();

    vector<std::string> Cuts;

       // Kst2Kpluspi0Resolved Cuts
    if(Kst2Kpluspi0Resolved){
        //don't copy old cone variables
        tree->SetBranchStatus("*_1.00",0);
        tree->SetBranchStatus("*_0.90",0);
        tree->SetBranchStatus("*_0.80",0);
        tree->SetBranchStatus("*_0.70",0);
        tree->SetBranchStatus("*_0.60",0);
        tree->SetBranchStatus("*_0.50",0);

        Cuts.push_back("K_plus_ProbNNk > 0.25");
        //Cuts.push_back("gamma1_PT_DTF < 4000");
        //Cuts.push_back("gamma2_PT_DTF < 4000");
        //Cuts.push_back("log(1.0-B_plus_DIRA_OWNPV) > -20");
    }
        // Kst2Kpluspi0Merged cuts
    if(Kst2Kpluspi0Merged){
        //don't copy old cone variables
        tree->SetBranchStatus("*_1.00",0);
        tree->SetBranchStatus("*_0.90",0);
        tree->SetBranchStatus("*_0.80",0);
        tree->SetBranchStatus("*_0.70",0);
        tree->SetBranchStatus("*_0.60",0);
        tree->SetBranchStatus("*_0.50",0);

        //#TODO

    }
        // Kst2Kspiplus cuts
    if(Kst2Kspiplus){
        //#TODO
    }



    std::string AllCuts = "";

    if (Cuts.size()==0){
        cout << "[ERROR]\t No cuts selected, nothing to do." << endl;
        return 0;
    }

    AllCuts = AllCuts + *Cuts.begin();
    cout << "[INFO]\t\t Appling cuts for year " << year << ":" << endl;
    cout << "\t\t\t" << AllCuts << endl;
    for (vector<std::string>::iterator t = Cuts.begin()+1;
                                       t != Cuts.end();
                                       ++t){
        AllCuts = AllCuts + " && " + *t ;
        cout << "\t\t\t" << *t << endl;
    }
    Cuts.clear();
    TTree* NewTree = tree->CopyTree( AllCuts.c_str() );

    output_new->cd();
    NewTree->Write();

    input->cd();

    NewTree->Write("",TObject::kOverwrite);

    if (useOld) output_old->Close();
    input->Close();
    output_new->Close();

    cout << "[INFO]\t\t New tree created." << endl;
    return 1;

}

Int_t sWeightSelectionAllData(Int_t Run = 1, bool useOld = false) {

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

}
Int_t sWeightSelectionAllMC(Int_t Run = 1, bool useOld = false) {

    for (auto& year : yearsMC(false, false, Run)){
        if ( sWeightSelection(useOld, year, "down", true, false, false) == 0) return 0;
        if ( sWeightSelection(useOld, year, "up",   true, false, false) == 0) return 0;
    }
    return 1;

}
Int_t sWeightSelectionAllPHSP(Int_t Run = 1, bool useOld = false) {

    for (auto& year : yearsMC(false, true, Run)){
        if ( sWeightSelection(useOld, year, "down", true, false, true) == 0) return 0;
        if ( sWeightSelection(useOld, year, "up",   true, false, true) == 0) return 0;
    }
    return 1;

}
Int_t sWeightSelectionAllRefChannel(Int_t Run = 1, bool useOld = false) {

    for (auto& year : yearsMC(true,false, Run)){
        if ( sWeightSelection(useOld, year, "down", true, true, false) == 0) return 0;
        if ( sWeightSelection(useOld, year, "up",   true, true, false) == 0) return 0;
    }
    return 1;

}
Int_t sWeightSelectionAll(Int_t Run = 1, bool useOld = false) {    

    if (sWeightSelectionAllData(Run,useOld)==0) return 0;
    if (sWeightSelectionAllMC(Run,useOld)==0) return 0;
    if (sWeightSelectionAllRefChannel(Run,useOld)==0) return 0;
    if (sWeightSelectionAllPHSP(Run,useOld)==0) return 0;
    return 1;

}

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

    if (PHSP) MC = true; //phase-space is always MC
    TChain* tree;

    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());
    }


   string output_path = GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,false);
   replace(output_path,".root","_StripForFit.root");
   // cout << output_path << endl;

   if(tree->GetEntries() == 0){
       //If it  complains while running, it's just trying to read also from David's preselected files
       std::cout << std::endl << "[ERROR]\t\tNo entries found in TTree from tfile " + GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample)+". Exit program!" << std::endl;
       std::cout << std::endl << "[ERROR]\t\tExit program!" << std::endl;
       return 0;
   }


   TFile *output = new TFile(output_path.c_str(),"RECREATE");

   //Get B mass
   Double_t B_plus_M;
   Double_t B_plus_M_DTF;
   Double_t B_plus_MERR_DTF;

   Float_t f_B_plus_M_DTF[100];
   Float_t f_B_plus_MERR_DTF[100];

   tree->SetBranchStatus("*",0);
   tree->SetBranchStatus("B_plus_M",1);
   tree->SetBranchStatus("B_plus_DTF_M",1);
   tree->SetBranchStatus("B_plus_DTF_MERR",1);

   tree -> SetBranchAddress( "B_plus_M" , &B_plus_M);
   tree -> SetBranchAddress( "B_plus_DTF_M" , &f_B_plus_M_DTF);
   tree -> SetBranchAddress( "B_plus_DTF_MERR" , &f_B_plus_MERR_DTF);

   //Get Q2
   tree->SetBranchStatus("mu_minus_PX",1);
   tree->SetBranchStatus("mu_minus_PY",1);
   tree->SetBranchStatus("mu_minus_PZ",1);
   tree->SetBranchStatus("mu_minus_PE",1);

   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_0_PX",1);
   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_0_PY",1);
   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_0_PZ",1);
   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_0_PE",1);

   tree->SetBranchStatus("mu_plus_PX",1);
   tree->SetBranchStatus("mu_plus_PY",1);
   tree->SetBranchStatus("mu_plus_PZ",1);
   tree->SetBranchStatus("mu_plus_PE",1);

   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_PX",1);
   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_PY",1);
   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_PZ",1);
   tree->SetBranchStatus("B_plus_DTF_J_psi_1S_muminus_PE",1);

   if (Kst2Kspiplus)tree->SetBranchStatus("B_plus_DTF_Kst_892_plus_piplus_ID",1);
   else tree->SetBranchStatus("B_plus_DTF_Kst_892_plus_Kplus_ID",1);


   Double_t mu_minus_PX     = 0.;
   Double_t mu_minus_PY     = 0.;
   Double_t mu_minus_PZ     = 0.;
   Double_t mu_minus_PE     = 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];

   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.;

      ///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.;

   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];

   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.;


  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( "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 );


   // Get Bplus charge to assign muplus and muminus correctly
   // Not really necessary, but why not
   //Not available for nonDTF at the moment!!!

   bool IsBplus = true;
   Double_t f_K_plus_ID_DTF[100];
   Double_t f_pi_plus_ID_DTF[100];

   if (Kst2Kspiplus) tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_ID" , &f_pi_plus_ID_DTF );
   else              tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_ID" , &f_K_plus_ID_DTF );


    //Define new branches
   TTree * NewTree = tree->CloneTree(0);
   NewTree->Branch("B_plus_M_DTF",        &B_plus_M_DTF,       "B_plus_M_DTF/D");
   NewTree->Branch("B_plus_MERR_DTF",     &B_plus_MERR_DTF,    "B_plus_MERR_DTF/D");

   unsigned int nDiMuonMassBin = 0;
   NewTree->Branch("nDiMuonMassBin", &nDiMuonMassBin, "nDiMuonMassBin/I");

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

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

   //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

   TLorentzVector LorVec_mu_minus;
   TLorentzVector LorVec_mu_plus;
   TLorentzVector LorVec_DiMuon;

       // Kst2Kpluspi0Resolved Cuts
    if(Kst2Kpluspi0Resolved){
        //don't copy old cone variables

        for (int evt = 0; evt < tree->GetEntries(); evt++){
            tree->GetEntry(evt);

            bool IsBplus = true;

            B_plus_M_DTF = (Double_t) f_B_plus_M_DTF[0];
            B_plus_MERR_DTF = (Double_t) f_B_plus_MERR_DTF[0];

            if (UseDTF && B_plus_M_DTF<cut_B_plus_M_low) continue;
            if (UseDTF && B_plus_M_DTF>cut_B_plus_M_high) continue;

            if (!UseDTF && B_plus_M<cut_B_plus_M_low) continue;
            if (!UseDTF && B_plus_M>cut_B_plus_M_high) continue;

            // --------- Get Q2 bin --------- //
            nDiMuonMassBin = 0;
            Is_Phi_JPsi_Psi2S = 0;
            if(Kst2Kspiplus) IsBplus = f_pi_plus_ID_DTF[0] > 0;
            else IsBplus = f_K_plus_ID_DTF[0] > 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_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];
            }
            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_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];
            }
            if (UseDTF){
                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);
                }
                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_minus_PE_DTF);
                }
            }
            else{
                if (UsePDGmIn4Dvectors){
                    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);
                }
                else{
                    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);
                }
            }

            LorVec_DiMuon = LorVec_mu_minus + LorVec_mu_plus;
            dQ2 = LorVec_DiMuon.M2();
            Is_Phi_JPsi_Psi2S = 0;
            //if(dQ2< f_BinBoundaries[0] || dQ2 > f_BinBoundaries[18])continue;
            while(dQ2 >= f_BinBoundaries[nDiMuonMassBin]){
                nDiMuonMassBin++;
            }
            if(nDiMuonMassBin == 2 || nDiMuonMassBin == 9 || nDiMuonMassBin == 12){
                Is_Phi_JPsi_Psi2S = 1;
            }
            ///////////////////////////////////////

            NewTree->Fill();
        }
    }
        // Kst2Kspiplus cuts
    if(Kst2Kspiplus){
        //#TODO
    }


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

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

}
int copyForFitAllData(){
    for (auto& year : yearsData(12)){
        if (copyForFit(year,"down",false,false,false)==0) return 0;
        if (copyForFit(year,"up",  false,false,false)==0) return 0;
    }
    return 1;

}
int copyForFitAllSig(){

    for (auto& year : yearsMC(false,false,12)){
        if (copyForFit(year,"down",true,false,false)==0) return 0;
        if (copyForFit(year,"up",  true,false,false)==0) return 0;
    }
    return 1;

}
int copyForFitAllRef(){

    for (auto& year : yearsMC(true,false,12)){
        if (copyForFit(year,"down",true,true,false)==0) return 0;
        if (copyForFit(year,"up",  true,true,false)==0) return 0;
    }
    return 1;

}
int copyForFitAllPHSP(){

    for (auto& year : yearsMC(false,true,12)){
        if (copyForFit(year,"down",true,false,false)==0) return 0;
        if (copyForFit(year,"up",  true,false,false)==0) return 0;
    }
    return 1;

}

int eventType(bool ReferenceChannel = false, bool PHSP = false){
    if (ReferenceChannel && PHSP) return 0;
    if(Kst2Kspiplus){
        if(ReferenceChannel) return 12145102;
        else if (PHSP) return 12115179;
        else return 12115102;
    }
    if(Kst2Kpluspi0Resolved){
        if(ReferenceChannel) return 12143401;
        else if (PHSP) return 12113446;
        else return 12113100;
    }
    return 0;
}
string eventTypeStr(bool ReferenceChannel = false, bool PHSP = false){
    if (ReferenceChannel && PHSP) return 0;
    if(Kst2Kspiplus){
        if(ReferenceChannel) return "12145102";
        else if (PHSP) return "12115179";
        else return "12115102";
    }
    if(Kst2Kpluspi0Resolved){
        if(ReferenceChannel) return "12143401";
        else if (PHSP) return "12113446";
        else return "12113100";
    }
    return 0;
}

int copyForPID(std::string year = "2011", std::string magnet = "down", bool ReferenceChannel = false, bool PHSP = false) {

    cout << "[INFO]\t\tCopy tuple for PID checks: " << year << magnet << endl;

    TChain* tree = new TChain();

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

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

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

    if(tree->GetEntries() > 0){
      std::cout << "[LOAD]\t\tSuccessfully loaded original TTree with " << tree->GetEntries() << " entries." << std::endl;
    }
    else
      return 0;

    string output_path = "";
    if(Kst2Kpluspi0Resolved)output_path = "/home/lhcb/kopecna/public/PID/Bu_Kstmumu_Kpi0_gammagamma_"+year+magnet+"_"+eventTypeStr(ReferenceChannel,PHSP)+".root";
    if(Kst2Kspiplus)        output_path = "/home/he/dgerick/Public/PID/Bu_Kstmumu_KSpi_pipi_"      +year+magnet+"_"+eventTypeStr(ReferenceChannel,PHSP)+".root";

    if(ReferenceChannel) replace(output_path,"mumu","Jpsi_mumu");


   TFile *output = new TFile(output_path.c_str(),"RECREATE");

   ///////////////////////////////////////
   //       Get particle variables      //
   ///////////////////////////////////////

   // ----- K+pi0 ----- //
   Double_t K_plus_P;       //Kst_0_K_0
   Double_t K_plus_PT;
   Int_t K_plus_TRUEID;
   Double_t K_plus_PIDe;
   Double_t K_plus_PIDmu;
   Double_t K_plus_PIDK;
   Double_t K_plus_PIDp;
   Double_t K_plus_ProbNNe;
   Double_t K_plus_ProbNNmu;
   Double_t K_plus_ProbNNk;
   Double_t K_plus_ProbNNp;
   Double_t K_plus_ProbNNpi;

   Double_t pi_zero_P;      //Kst_0_Pi0_0
   Double_t pi_zero_PT;
   Int_t pi_zero_TRUEID;
   Double_t pi_zero_M;
   Double_t pi_zero_isPhoton;
   Double_t pi_zero_isNotH;
   Double_t pi_zero_isNotE;


   Double_t Gamma_0_P;      //Pi0_0_Gamma_0
   Double_t Gamma_0_PT;
   Int_t Gamma_0_TRUEID;
   Double_t Gamma_0_isPhoton;
   Double_t Gamma_0_isNotH;
   Double_t Gamma_0_isNotE;

   Double_t Gamma_1_P;      //Pi0_0_Gamma_1
   Double_t Gamma_1_PT;
   Int_t Gamma_1_TRUEID;
   Double_t Gamma_1_isPhoton;
   Double_t Gamma_1_isNotH;
   Double_t Gamma_1_isNotE;


   // ----- Kspi+ ----- //

   Double_t K_short_P;      //Kst_0_K_0
   Double_t K_short_PT;
   Int_t K_short_TRUEID;

   Double_t pi_plus_P;      //Kst_0_Pi_0
   Double_t pi_plus_PT;
   Int_t pi_plus_TRUEID;
   Double_t pi_plus_PIDe;
   Double_t pi_plus_PIDmu;
   Double_t pi_plus_PIDK;
   Double_t pi_plus_PIDp;
   Double_t pi_plus_ProbNNe;
   Double_t pi_plus_ProbNNmu;
   Double_t pi_plus_ProbNNk;
   Double_t pi_plus_ProbNNp;
   Double_t pi_plus_ProbNNpi;

   Double_t Ks_pi_plus_P;   //Ks_0_Pi_0
   Double_t Ks_pi_plus_PT;
   Int_t Ks_pi_plus_TRUEID;
   Double_t Ks_pi_plus_PIDe;
   Double_t Ks_pi_plus_PIDmu;
   Double_t Ks_pi_plus_PIDK;
   Double_t Ks_pi_plus_PIDp;
   Double_t Ks_pi_plus_ProbNNe;
   Double_t Ks_pi_plus_ProbNNmu;
   Double_t Ks_pi_plus_ProbNNk;
   Double_t Ks_pi_plus_ProbNNp;
   Double_t Ks_pi_plus_ProbNNpi;

   Double_t Ks_pi_minus_P;  //Ks_0_Pi_1
   Double_t Ks_pi_minus_PT;
   Int_t Ks_pi_minus_TRUEID;
   Double_t Ks_pi_minus_PIDe;
   Double_t Ks_pi_minus_PIDmu;
   Double_t Ks_pi_minus_PIDK;
   Double_t Ks_pi_minus_PIDp;
   Double_t Ks_pi_minus_ProbNNe;
   Double_t Ks_pi_minus_ProbNNmu;
   Double_t Ks_pi_minus_ProbNNk;
   Double_t Ks_pi_minus_ProbNNp;
   Double_t Ks_pi_minus_ProbNNpi;

   // ---- mumu ----- //

   Double_t mu_plus_P;   //B_0_Mu_0
   Double_t mu_plus_PT;
   Int_t mu_plus_TRUEID;
   Double_t mu_plus_PIDe;
   Double_t mu_plus_PIDmu;
   Double_t mu_plus_PIDK;
   Double_t mu_plus_PIDp;
   Double_t mu_plus_ProbNNe;
   Double_t mu_plus_ProbNNmu;
   Double_t mu_plus_ProbNNk;
   Double_t mu_plus_ProbNNp;
   Double_t mu_plus_ProbNNpi;

   Double_t mu_minus_P;  //B_0_Mu_1
   Double_t mu_minus_PT;
   Int_t mu_minus_TRUEID;
   Double_t mu_minus_PIDe;
   Double_t mu_minus_PIDmu;
   Double_t mu_minus_PIDK;
   Double_t mu_minus_PIDp;
   Double_t mu_minus_ProbNNe;
   Double_t mu_minus_ProbNNmu;
   Double_t mu_minus_ProbNNk;
   Double_t mu_minus_ProbNNp;
   Double_t mu_minus_ProbNNpi;


   ///////////////////////////////////////
   //            Load branches          //
   ///////////////////////////////////////

   tree->SetBranchStatus("*",0);

   // ----- K+pi0 ----- //
   if(Kst2Kpluspi0Resolved){

       tree->SetBranchStatus("K_plus_P",1);       //Kst_0_K_0
       tree->SetBranchStatus("K_plus_PT",1);
       tree->SetBranchStatus("K_plus_TRUEID",1);
       tree->SetBranchStatus("K_plus_PIDe",1);
       tree->SetBranchStatus("K_plus_PIDmu",1);
       tree->SetBranchStatus("K_plus_PIDK",1);
       tree->SetBranchStatus("K_plus_PIDp",1);
       tree->SetBranchStatus("K_plus_ProbNNe",1);
       tree->SetBranchStatus("K_plus_ProbNNmu",1);
       tree->SetBranchStatus("K_plus_ProbNNk",1);
       tree->SetBranchStatus("K_plus_ProbNNp",1);
       tree->SetBranchStatus("K_plus_ProbNNpi",1);

       tree->SetBranchAddress("K_plus_P",&K_plus_P);      //Kst_0_K_0
       tree->SetBranchAddress("K_plus_PT",&K_plus_PT);
       tree->SetBranchAddress("K_plus_TRUEID",&K_plus_TRUEID);
       tree->SetBranchAddress("K_plus_PIDe",&K_plus_PIDe);
       tree->SetBranchAddress("K_plus_PIDmu",&K_plus_PIDmu);
       tree->SetBranchAddress("K_plus_PIDK",&K_plus_PIDK);
       tree->SetBranchAddress("K_plus_PIDp",&K_plus_PIDp);
       tree->SetBranchAddress("K_plus_ProbNNe",&K_plus_ProbNNe);
       tree->SetBranchAddress("K_plus_ProbNNmu",&K_plus_ProbNNmu);
       tree->SetBranchAddress("K_plus_ProbNNk",&K_plus_ProbNNk);
       tree->SetBranchAddress("K_plus_ProbNNp",&K_plus_ProbNNp);
       tree->SetBranchAddress("K_plus_ProbNNpi",&K_plus_ProbNNpi);

       tree->SetBranchStatus("pi_zero_resolved_P",1);      //Kst_0_Pi0_0
       tree->SetBranchStatus("pi_zero_resolved_PT",1);
       tree->SetBranchStatus("pi_zero_resolved_TRUEID",1);
       tree->SetBranchStatus("pi_zero_resolved_M",1);
       //tree->SetBranchStatus("pi_zero_resolved_isPhoton",1);
       //tree->SetBranchStatus("pi_zero_resolved_isNotH",1);
       //tree->SetBranchStatus("pi_zero_resolved_isNotE",1);

       tree->SetBranchAddress("pi_zero_resolved_P",&pi_zero_P);     //Kst_0_Pi0_0
       tree->SetBranchAddress("pi_zero_resolved_PT",&pi_zero_PT);
       tree->SetBranchAddress("pi_zero_resolved_TRUEID",&pi_zero_TRUEID);
       tree->SetBranchAddress("pi_zero_resolved_M",&pi_zero_M);
       //tree->SetBranchAddress("pi_zero_resolved_isPhoton",&pi_zero_isPhoton);
       //tree->SetBranchAddress("pi_zero_resolved_isNotH",&pi_zero_isNotH);
       //tree->SetBranchAddress("pi_zero_resolved_isNotE",&pi_zero_isNotE);

       tree->SetBranchStatus("gamma1_P",1);      //Pi0_0_Gamma_0
       tree->SetBranchStatus("gamma1_PT",1);
       tree->SetBranchStatus("gamma1_TRUEID",1);
       //tree->SetBranchStatus("gamma1_isPhoton",1);
       //tree->SetBranchStatus("gamma1_isNotH",1);
       //tree->SetBranchStatus("gamma1_isNotE",1);

       tree->SetBranchStatus("gamma2_P",1);      //Pi0_0_Gamma_1
       tree->SetBranchStatus("gamma2_PT",1);
       tree->SetBranchStatus("gamma2_TRUEID",1);
       //tree->SetBranchStatus("gamma2_isPhoton",1);
       //tree->SetBranchStatus("gamma2_isNotH",1);
       //tree->SetBranchStatus("gamma2_isNotE",1);

       tree->SetBranchAddress("gamma1_P",&Gamma_0_P);      //Pi0_0_Gamma_0
       tree->SetBranchAddress("gamma1_PT",&Gamma_0_PT);
       tree->SetBranchAddress("gamma1_TRUEID",&Gamma_0_TRUEID);
       //tree->SetBranchAddress("gamma1_isPhoton",&Gamma_0_isPhoton);
       //tree->SetBranchAddress("gamma1_isNotH",&Gamma_0_isNotH);
       //tree->SetBranchAddress("gamma1_isNotE",&Gamma_0_isNotE);

       tree->SetBranchAddress("gamma2_P",&Gamma_1_P);      //Pi0_0_Gamma_1
       tree->SetBranchAddress("gamma2_PT",&Gamma_1_PT);
       tree->SetBranchAddress("gamma2_TRUEID",&Gamma_1_TRUEID);
       //tree->SetBranchAddress("gamma2_isPhoton",&Gamma_1_isPhoton);
       //tree->SetBranchAddress("gamma2_isNotH",&Gamma_1_isNotH);
       //tree->SetBranchAddress("gamma2_isNotE",&Gamma_1_isNotE);

   }

   // ----- Kspi+ ----- //
   if(Kst2Kspiplus){

       tree->SetBranchStatus("K_short_P",1);      //Kst_0_K_0
       tree->SetBranchStatus("K_short_PT",1);
       tree->SetBranchStatus("K_short_TRUEID",1);

       tree->SetBranchStatus("pi_plus_P",1);      //Kst_0_Pi_0
       tree->SetBranchStatus("pi_plus_PT",1);
       tree->SetBranchStatus("pi_plus_TRUEID",1);
       tree->SetBranchStatus("pi_plus_PIDe",1);
       tree->SetBranchStatus("pi_plus_PIDmu",1);
       tree->SetBranchStatus("pi_plus_PIDK",1);
       tree->SetBranchStatus("pi_plus_PIDp",1);
       tree->SetBranchStatus("pi_plus_ProbNNe",1);
       tree->SetBranchStatus("pi_plus_ProbNNmu",1);
       tree->SetBranchStatus("pi_plus_ProbNNk",1);
       tree->SetBranchStatus("pi_plus_ProbNNp",1);
       tree->SetBranchStatus("pi_plus_ProbNNpi",1);

       tree->SetBranchStatus("Ks_pi_plus_P",1);   //Ks_0_Pi_0
       tree->SetBranchStatus("Ks_pi_plus_PT",1);
       tree->SetBranchStatus("Ks_pi_plus_TRUEID",1);
       tree->SetBranchStatus("Ks_pi_plus_PIDe",1);
       tree->SetBranchStatus("Ks_pi_plus_PIDmu",1);
       tree->SetBranchStatus("Ks_pi_plus_PIDK",1);
       tree->SetBranchStatus("Ks_pi_plus_PIDp",1);
       tree->SetBranchStatus("Ks_pi_plus_ProbNNe",1);
       tree->SetBranchStatus("Ks_pi_plus_ProbNNmu",1);
       tree->SetBranchStatus("Ks_pi_plus_ProbNNk",1);
       tree->SetBranchStatus("Ks_pi_plus_ProbNNp",1);
       tree->SetBranchStatus("Ks_pi_plus_ProbNNpi",1);

       tree->SetBranchStatus("Ks_pi_minus_P",1);  //Ks_0_Pi_1
       tree->SetBranchStatus("Ks_pi_minus_PT",1);
       tree->SetBranchStatus("Ks_pi_minus_TRUEID",1);
       tree->SetBranchStatus("Ks_pi_minus_PIDe",1);
       tree->SetBranchStatus("Ks_pi_minus_PIDmu",1);
       tree->SetBranchStatus("Ks_pi_minus_PIDK",1);
       tree->SetBranchStatus("Ks_pi_minus_PIDp",1);
       tree->SetBranchStatus("Ks_pi_minus_ProbNNe",1);
       tree->SetBranchStatus("Ks_pi_minus_ProbNNmu",1);
       tree->SetBranchStatus("Ks_pi_minus_ProbNNk",1);
       tree->SetBranchStatus("Ks_pi_minus_ProbNNp",1);
       tree->SetBranchStatus("Ks_pi_minus_ProbNNpi",1);

       tree->SetBranchAddress("K_short_P",&K_short_P);      //Kst_0_K_0
       tree->SetBranchAddress("K_short_PT",&K_short_PT);
       tree->SetBranchAddress("K_short_TRUEID",&K_short_TRUEID);;

       tree->SetBranchAddress("pi_plus_P",&pi_plus_P);      //Kst_0_Pi_0
       tree->SetBranchAddress("pi_plus_PT",&pi_plus_PT);
       tree->SetBranchAddress("pi_plus_TRUEID",&pi_plus_TRUEID);
       tree->SetBranchAddress("pi_plus_PIDe",&pi_plus_PIDe);
       tree->SetBranchAddress("pi_plus_PIDmu",&pi_plus_PIDmu);
       tree->SetBranchAddress("pi_plus_PIDK",&pi_plus_PIDK);
       tree->SetBranchAddress("pi_plus_PIDp",&pi_plus_PIDp);
       tree->SetBranchAddress("pi_plus_ProbNNe",&pi_plus_ProbNNe);
       tree->SetBranchAddress("pi_plus_ProbNNmu",&pi_plus_ProbNNmu);
       tree->SetBranchAddress("pi_plus_ProbNNk",&pi_plus_ProbNNk);
       tree->SetBranchAddress("pi_plus_ProbNNp",&pi_plus_ProbNNp);
       tree->SetBranchAddress("pi_plus_ProbNNpi",&pi_plus_ProbNNpi);

       tree->SetBranchAddress("Ks_pi_plus_P",&Ks_pi_plus_P);   //Ks_0_Pi_0
       tree->SetBranchAddress("Ks_pi_plus_PT",&Ks_pi_plus_PT);
       tree->SetBranchAddress("Ks_pi_plus_TRUEID",&Ks_pi_plus_TRUEID);
       tree->SetBranchAddress("Ks_pi_plus_PIDe",&Ks_pi_plus_PIDe);
       tree->SetBranchAddress("Ks_pi_plus_PIDmu",&Ks_pi_plus_PIDmu);
       tree->SetBranchAddress("Ks_pi_plus_PIDK",&Ks_pi_plus_PIDK);
       tree->SetBranchAddress("Ks_pi_plus_PIDp",&Ks_pi_plus_PIDp);
       tree->SetBranchAddress("Ks_pi_plus_ProbNNe",&Ks_pi_plus_ProbNNe);
       tree->SetBranchAddress("Ks_pi_plus_ProbNNmu",&Ks_pi_plus_ProbNNmu);
       tree->SetBranchAddress("Ks_pi_plus_ProbNNk",&Ks_pi_plus_ProbNNk);
       tree->SetBranchAddress("Ks_pi_plus_ProbNNp",&Ks_pi_plus_ProbNNp);
       tree->SetBranchAddress("Ks_pi_plus_ProbNNpi",&Ks_pi_plus_ProbNNpi);

       tree->SetBranchAddress("Ks_pi_minus_P",&Ks_pi_minus_P);  //Ks_0_Pi_1
       tree->SetBranchAddress("Ks_pi_minus_PT",&Ks_pi_minus_PT);
       tree->SetBranchAddress("Ks_pi_minus_TRUEID",&Ks_pi_minus_TRUEID);
       tree->SetBranchAddress("Ks_pi_minus_PIDe",&Ks_pi_minus_PIDe);
       tree->SetBranchAddress("Ks_pi_minus_PIDmu",&Ks_pi_minus_PIDmu);
       tree->SetBranchAddress("Ks_pi_minus_PIDK",&Ks_pi_minus_PIDK);
       tree->SetBranchAddress("Ks_pi_minus_PIDp",&Ks_pi_minus_PIDp);
       tree->SetBranchAddress("Ks_pi_minus_ProbNNe",&Ks_pi_minus_ProbNNe);
       tree->SetBranchAddress("Ks_pi_minus_ProbNNmu",&Ks_pi_minus_ProbNNmu);
       tree->SetBranchAddress("Ks_pi_minus_ProbNNk",&Ks_pi_minus_ProbNNk);
       tree->SetBranchAddress("Ks_pi_minus_ProbNNp",&Ks_pi_minus_ProbNNp);
       tree->SetBranchAddress("Ks_pi_minus_ProbNNpi",&Ks_pi_minus_ProbNNpi);
   }

   // ---- mumu ----- //
   tree->SetBranchStatus("mu_plus_P",1);   //B_0_Mu_0
   tree->SetBranchStatus("mu_plus_PT",1);
   tree->SetBranchStatus("mu_plus_TRUEID",1);
   tree->SetBranchStatus("mu_plus_PIDe",1);
   tree->SetBranchStatus("mu_plus_PIDmu",1);
   tree->SetBranchStatus("mu_plus_PIDK",1);
   tree->SetBranchStatus("mu_plus_PIDp",1);
   tree->SetBranchStatus("mu_plus_ProbNNe",1);
   tree->SetBranchStatus("mu_plus_ProbNNmu",1);
   tree->SetBranchStatus("mu_plus_ProbNNk",1);
   tree->SetBranchStatus("mu_plus_ProbNNp",1);
   tree->SetBranchStatus("mu_plus_ProbNNpi",1);

   tree->SetBranchStatus("mu_minus_P",1);  //B_0_Mu_1
   tree->SetBranchStatus("mu_minus_PT",1);
   tree->SetBranchStatus("mu_minus_TRUEID",1);
   tree->SetBranchStatus("mu_minus_PIDe",1);
   tree->SetBranchStatus("mu_minus_PIDmu",1);
   tree->SetBranchStatus("mu_minus_PIDK",1);
   tree->SetBranchStatus("mu_minus_PIDp",1);
   tree->SetBranchStatus("mu_minus_ProbNNe",1);
   tree->SetBranchStatus("mu_minus_ProbNNmu",1);
   tree->SetBranchStatus("mu_minus_ProbNNk",1);
   tree->SetBranchStatus("mu_minus_ProbNNp",1);
   tree->SetBranchStatus("mu_minus_ProbNNpi",1);

   tree->SetBranchAddress("mu_plus_P",&mu_plus_P);   //B_0_Mu_0
   tree->SetBranchAddress("mu_plus_PT",&mu_plus_PT);
   tree->SetBranchAddress("mu_plus_TRUEID",&mu_plus_TRUEID);
   tree->SetBranchAddress("mu_plus_PIDe",&mu_plus_PIDe);
   tree->SetBranchAddress("mu_plus_PIDmu",&mu_plus_PIDmu);
   tree->SetBranchAddress("mu_plus_PIDK",&mu_plus_PIDK);
   tree->SetBranchAddress("mu_plus_PIDp",&mu_plus_PIDp);
   tree->SetBranchAddress("mu_plus_ProbNNe",&mu_plus_ProbNNe);
   tree->SetBranchAddress("mu_plus_ProbNNmu",&mu_plus_ProbNNmu);
   tree->SetBranchAddress("mu_plus_ProbNNk",&mu_plus_ProbNNk);
   tree->SetBranchAddress("mu_plus_ProbNNp",&mu_plus_ProbNNp);
   tree->SetBranchAddress("mu_plus_ProbNNpi",&mu_plus_ProbNNpi);

   tree->SetBranchAddress("mu_minus_P",&mu_minus_P);  //B_0_Mu_1
   tree->SetBranchAddress("mu_minus_PT",&mu_minus_PT);
   tree->SetBranchAddress("mu_minus_TRUEID",&mu_minus_TRUEID);
   tree->SetBranchAddress("mu_minus_PIDe",&mu_minus_PIDe);
   tree->SetBranchAddress("mu_minus_PIDmu",&mu_minus_PIDmu);
   tree->SetBranchAddress("mu_minus_PIDK",&mu_minus_PIDK);
   tree->SetBranchAddress("mu_minus_PIDp",&mu_minus_PIDp);
   tree->SetBranchAddress("mu_minus_ProbNNe",&mu_minus_ProbNNe);
   tree->SetBranchAddress("mu_minus_ProbNNmu",&mu_minus_ProbNNmu);
   tree->SetBranchAddress("mu_minus_ProbNNk",&mu_minus_ProbNNk);
   tree->SetBranchAddress("mu_minus_ProbNNp",&mu_minus_ProbNNp);
   tree->SetBranchAddress("mu_minus_ProbNNpi",&mu_minus_ProbNNpi);


   ///////////////////////////////////////
   //        Define new branches        //
   ///////////////////////////////////////

   TTree * NewTree = new TTree("DecayTree", "DecayTree");

   // ----- K+pi0 ----- //
   if(Kst2Kpluspi0Resolved){

       NewTree->Branch("Kst_0_K_0_P",&K_plus_P,"Kst_0_K_0_P/D");      //Kst_0_K_0
       NewTree->Branch("Kst_0_K_0_PT",&K_plus_PT,"Kst_0_K_0_PT/D");
       NewTree->Branch("Kst_0_K_0_TRUEID",&K_plus_TRUEID,"Kst_0_K_0_TRUEID/I");
       NewTree->Branch("Kst_0_K_0_PIDe",&K_plus_PIDe,"Kst_0_K_0_PIDe/D");
       NewTree->Branch("Kst_0_K_0_PIDmu",&K_plus_PIDmu,"Kst_0_K_0_PIDmu/D");
       NewTree->Branch("Kst_0_K_0_PIDK",&K_plus_PIDK,"Kst_0_K_0_PIDK/D");
       NewTree->Branch("Kst_0_K_0_PIDp",&K_plus_PIDp,"Kst_0_K_0_PIDp/D");
       NewTree->Branch("Kst_0_K_0_ProbNNe",&K_plus_ProbNNe,"Kst_0_K_0_ProbNNe/D");
       NewTree->Branch("Kst_0_K_0_ProbNNmu",&K_plus_ProbNNmu,"Kst_0_K_0_ProbNNmu/D");
       NewTree->Branch("Kst_0_K_0_ProbNNk",&K_plus_ProbNNk,"Kst_0_K_0_ProbNNk/D");
       NewTree->Branch("Kst_0_K_0_ProbNNp",&K_plus_ProbNNp,"Kst_0_K_0_ProbNNp/D");
       NewTree->Branch("Kst_0_K_0_ProbNNpi",&K_plus_ProbNNpi,"Kst_0_K_0_ProbNNpi/D");

       NewTree->Branch("Kst_0_Pi0_0_P",&pi_zero_P,"Kst_0_Pi0_0_P/D");     //Kst_0_Pi0_0
       NewTree->Branch("Kst_0_Pi0_0_PT",&pi_zero_PT,"Kst_0_Pi0_0_PT/D");
       NewTree->Branch("Kst_0_Pi0_0_TRUEID",&pi_zero_TRUEID,"Kst_0_Pi0_0_TRUEID/I");
       NewTree->Branch("Kst_0_Pi0_0_M",&pi_zero_M,"Kst_0_Pi0_0_M/D");
       //NewTree->Branch("Kst_0_Pi0_0_isPhoton",&pi_zero_isPhoton,"Kst_0_Pi0_0_isPhoton/D");
       //NewTree->Branch("Kst_0_Pi0_0_isNotH",&pi_zero_isNotH,"Kst_0_Pi0_0_isNotH/D");
       //NewTree->Branch("Kst_0_Pi0_0_isNotE",&pi_zero_isNotE,"Kst_0_Pi0_0_isNotE/D");

       NewTree->Branch("Pi0_0_Gamma_0_P",&Gamma_0_P,"Pi0_0_Gamma_0_P/D");      //Pi0_0_Gamma_0
       NewTree->Branch("Pi0_0_Gamma_0_PT",&Gamma_0_PT,"Pi0_0_Gamma_0_PT/D");
       NewTree->Branch("Pi0_0_Gamma_0_TRUEID",&Gamma_0_TRUEID,"Pi0_0_Gamma_0_TRUEID/I");
       //NewTree->Branch("Pi0_0_Gamma_0_isPhoton",&Gamma_0_isPhoton,"Pi0_0_Gamma_0_isPhoton/D");
       //NewTree->Branch("Pi0_0_Gamma_0_isNotH",&Gamma_0_isNotH,"Pi0_0_Gamma_0_isNotH/D");
       //NewTree->Branch("Pi0_0_Gamma_0_isNotE",&Gamma_0_isNotE,"Pi0_0_Gamma_0_isNotE/D");

       NewTree->Branch("Pi0_0_Gamma_1_P",&Gamma_1_P,"Pi0_0_Gamma_1_P/D");      //Pi0_0_Gamma_1
       NewTree->Branch("Pi0_0_Gamma_1_PT",&Gamma_1_PT,"Pi0_0_Gamma_1_PT/D");
       NewTree->Branch("Pi0_0_Gamma_1_TRUEID",&Gamma_1_TRUEID,"Pi0_0_Gamma_1_TRUEID/I");
       //NewTree->Branch("Pi0_0_Gamma_1_isPhoton",&Gamma_1_isPhoton,"Pi0_0_Gamma_1_isPhoton/D");
       //NewTree->Branch("Pi0_0_Gamma_1_isNotH",&Gamma_1_isNotH,"Pi0_0_Gamma_1_isNotH/D");
       //NewTree->Branch("Pi0_0_Gamma_1_isNotE",&Gamma_1_isNotE,"Pi0_0_Gamma_1_isNotE/D");

   }

   // ----- Kspi+ ----- //
   if(Kst2Kspiplus){
       NewTree->Branch("Kst_0_K_0_P",&K_short_P,"Kst_0_K_0_P/D");      //Kst_0_K_0
       NewTree->Branch("Kst_0_K_0_PT",&K_short_PT,"Kst_0_K_0_PT/D");
       NewTree->Branch("Kst_0_K_0_TRUEID",&K_short_TRUEID,"Kst_0_K_0_TRUEID/I");

       NewTree->Branch("Kst_0_Pi_0_P",&pi_plus_P,"Kst_0_Pi_0_P/D");      //Kst_0_Pi_0
       NewTree->Branch("Kst_0_Pi_0_PT",&pi_plus_PT,"Kst_0_Pi_0_PT/D");
       NewTree->Branch("Kst_0_Pi_0_TRUEID",&pi_plus_TRUEID,"Kst_0_Pi_0_TRUEID/I");
       NewTree->Branch("Kst_0_Pi_0_PIDe",&pi_plus_PIDe,"Kst_0_Pi_0_PIDe/D");
       NewTree->Branch("Kst_0_Pi_0_PIDmu",&pi_plus_PIDmu,"Kst_0_Pi_0_PIDmu/D");
       NewTree->Branch("Kst_0_Pi_0_PIDK",&pi_plus_PIDK,"Kst_0_Pi_0_PIDK/D");
       NewTree->Branch("Kst_0_Pi_0_PIDp",&pi_plus_PIDp,"Kst_0_Pi_0_PIDp/D");
       NewTree->Branch("Kst_0_Pi_0_ProbNNe",&pi_plus_ProbNNe,"Kst_0_Pi_0_ProbNNe/D");
       NewTree->Branch("Kst_0_Pi_0_ProbNNmu",&pi_plus_ProbNNmu,"Kst_0_Pi_0_ProbNNmu/D");
       NewTree->Branch("Kst_0_Pi_0_ProbNNk",&pi_plus_ProbNNk,"Kst_0_Pi_0_ProbNNk/D");
       NewTree->Branch("Kst_0_Pi_0_ProbNNp",&pi_plus_ProbNNp,"Kst_0_Pi_0_ProbNNp/D");
       NewTree->Branch("Kst_0_Pi_0_ProbNNpi",&pi_plus_ProbNNpi,"Kst_0_Pi_0_ProbNNpi/D");

       NewTree->Branch("K_0_Pi_0_P",&Ks_pi_plus_P,"K_0_Pi_0_P/D");   //Ks_0_Pi_0
       NewTree->Branch("K_0_Pi_0_PT",&Ks_pi_plus_PT,"K_0_Pi_0_PT/D");
       NewTree->Branch("K_0_Pi_0_TRUEID",&Ks_pi_plus_TRUEID,"K_0_Pi_0_TRUEID/I");
       NewTree->Branch("K_0_Pi_0_PIDe",&Ks_pi_plus_PIDe,"K_0_Pi_0_PIDe/D");
       NewTree->Branch("K_0_Pi_0_PIDmu",&Ks_pi_plus_PIDmu,"K_0_Pi_0_PIDmu/D");
       NewTree->Branch("K_0_Pi_0_PIDK",&Ks_pi_plus_PIDK,"K_0_Pi_0_PIDK/D");
       NewTree->Branch("K_0_Pi_0_PIDp",&Ks_pi_plus_PIDp,"K_0_Pi_0_PIDp/D");
       NewTree->Branch("K_0_Pi_0_ProbNNe",&Ks_pi_plus_ProbNNe,"K_0_Pi_0_ProbNNe/D");
       NewTree->Branch("K_0_Pi_0_ProbNNmu",&Ks_pi_plus_ProbNNmu,"K_0_Pi_0_ProbNNmu/D");
       NewTree->Branch("K_0_Pi_0_ProbNNk",&Ks_pi_plus_ProbNNk,"K_0_Pi_0_ProbNNk/D");
       NewTree->Branch("K_0_Pi_0_ProbNNp",&Ks_pi_plus_ProbNNp,"K_0_Pi_0_ProbNNp/D");
       NewTree->Branch("K_0_Pi_0_ProbNNpi",&Ks_pi_plus_ProbNNpi,"K_0_Pi_0_ProbNNpi/D");

       NewTree->Branch("K_0_Pi_1_P",&Ks_pi_minus_P,"K_0_Pi_1_P/D");  //Ks_0_Pi_1
       NewTree->Branch("K_0_Pi_1_PT",&Ks_pi_minus_PT,"K_0_Pi_1_PT/D");
       NewTree->Branch("K_0_Pi_1_TRUEID",&Ks_pi_minus_TRUEID,"K_0_Pi_1_TRUEID/I");
       NewTree->Branch("K_0_Pi_1_PIDe",&Ks_pi_minus_PIDe,"K_0_Pi_1_PIDe/D");
       NewTree->Branch("K_0_Pi_1_PIDmu",&Ks_pi_minus_PIDmu,"K_0_Pi_1_PIDmu/D");
       NewTree->Branch("K_0_Pi_1_PIDK",&Ks_pi_minus_PIDK,"K_0_Pi_1_PIDK/D");
       NewTree->Branch("K_0_Pi_1_PIDp",&Ks_pi_minus_PIDp,"K_0_Pi_1_PIDp/D");
       NewTree->Branch("K_0_Pi_1_ProbNNe",&Ks_pi_minus_ProbNNe,"K_0_Pi_1_ProbNNe/D");
       NewTree->Branch("K_0_Pi_1_ProbNNmu",&Ks_pi_minus_ProbNNmu,"K_0_Pi_1_ProbNNmu/D");
       NewTree->Branch("K_0_Pi_1_ProbNNk",&Ks_pi_minus_ProbNNk,"K_0_Pi_1_ProbNNk/D");
       NewTree->Branch("K_0_Pi_1_ProbNNp",&Ks_pi_minus_ProbNNp,"K_0_Pi_1_ProbNNp/D");
       NewTree->Branch("K_0_Pi_1_ProbNNpi",&Ks_pi_minus_ProbNNpi,"K_0_Pi_1_ProbNNpi/D");
   }

   // ---- mumu ----- //
   if (ReferenceChannel){
       NewTree->Branch("Jpsi_0_Mu_0_P",&mu_plus_P,"Jpsi_0_Mu_0_P/D");   //Jpsi_0_Mu_0
       NewTree->Branch("Jpsi_0_Mu_0_PT",&mu_plus_PT,"Jpsi_0_Mu_0_PT/D");
       NewTree->Branch("Jpsi_0_Mu_0_TRUEID",&mu_plus_TRUEID,"Jpsi_0_Mu_0_TRUEID/I");
       NewTree->Branch("Jpsi_0_Mu_0_PIDe",&mu_plus_PIDe,"Jpsi_0_Mu_0_PIDe/D");
       NewTree->Branch("Jpsi_0_Mu_0_PIDmu",&mu_plus_PIDmu,"Jpsi_0_Mu_0_PIDmu/D");
       NewTree->Branch("Jpsi_0_Mu_0_PIDK",&mu_plus_PIDK,"Jpsi_0_Mu_0_PIDK/D");
       NewTree->Branch("Jpsi_0_Mu_0_PIDp",&mu_plus_PIDp,"Jpsi_0_Mu_0_PIDp/D");
       NewTree->Branch("Jpsi_0_Mu_0_ProbNNe",&mu_plus_ProbNNe,"Jpsi_0_Mu_0_ProbNNe/D");
       NewTree->Branch("Jpsi_0_Mu_0_ProbNNmu",&mu_plus_ProbNNmu,"Jpsi_0_Mu_0_ProbNNmu/D");
       NewTree->Branch("Jpsi_0_Mu_0_ProbNNk",&mu_plus_ProbNNk,"Jpsi_0_Mu_0_ProbNNk/D");
       NewTree->Branch("Jpsi_0_Mu_0_ProbNNp",&mu_plus_ProbNNp,"Jpsi_0_Mu_0_ProbNNp/D");
       NewTree->Branch("Jpsi_0_Mu_0_ProbNNpi",&mu_plus_ProbNNpi,"Jpsi_0_Mu_0_ProbNNpi/D");

       NewTree->Branch("Jpsi_0_Mu_1_P",&mu_minus_P,"Jpsi_0_Mu_1_P/D");  //Jpsi_0_Mu_1
       NewTree->Branch("Jpsi_0_Mu_1_PT",&mu_minus_PT,"Jpsi_0_Mu_1_PT/D");
       NewTree->Branch("Jpsi_0_Mu_1_TRUEID",&mu_minus_TRUEID,"Jpsi_0_Mu_1_TRUEID/I");
       NewTree->Branch("Jpsi_0_Mu_1_PIDe",&mu_minus_PIDe,"Jpsi_0_Mu_1_PIDe/D");
       NewTree->Branch("Jpsi_0_Mu_1_PIDmu",&mu_minus_PIDmu,"Jpsi_0_Mu_1_PIDmu/D");
       NewTree->Branch("Jpsi_0_Mu_1_PIDK",&mu_minus_PIDK,"Jpsi_0_Mu_1_PIDK/D");
       NewTree->Branch("Jpsi_0_Mu_1_PIDp",&mu_minus_PIDp,"Jpsi_0_Mu_1_PIDp/D");
       NewTree->Branch("Jpsi_0_Mu_1_ProbNNe",&mu_minus_ProbNNe,"Jpsi_0_Mu_1_ProbNNe/D");
       NewTree->Branch("Jpsi_0_Mu_1_ProbNNmu",&mu_minus_ProbNNmu,"Jpsi_0_Mu_1_ProbNNmu/D");
       NewTree->Branch("Jpsi_0_Mu_1_ProbNNk",&mu_minus_ProbNNk,"Jpsi_0_Mu_1_ProbNNk/D");
       NewTree->Branch("Jpsi_0_Mu_1_ProbNNp",&mu_minus_ProbNNp,"Jpsi_0_Mu_1_ProbNNp/D");
       NewTree->Branch("Jpsi_0_Mu_1_ProbNNpi",&mu_minus_ProbNNpi,"Jpsi_0_Mu_1_ProbNNpi/D");
   }
   else{
       NewTree->Branch("B_0_Mu_0_P",&mu_plus_P,"B_0_Mu_0_P/D");   //B_0_Mu_0
       NewTree->Branch("B_0_Mu_0_PT",&mu_plus_PT,"B_0_Mu_0_PT/D");
       NewTree->Branch("B_0_Mu_0_TRUEID",&mu_plus_TRUEID,"B_0_Mu_0_TRUEID/I");
       NewTree->Branch("B_0_Mu_0_PIDe",&mu_plus_PIDe,"B_0_Mu_0_PIDe/D");
       NewTree->Branch("B_0_Mu_0_PIDmu",&mu_plus_PIDmu,"B_0_Mu_0_PIDmu/D");
       NewTree->Branch("B_0_Mu_0_PIDK",&mu_plus_PIDK,"B_0_Mu_0_PIDK/D");
       NewTree->Branch("B_0_Mu_0_PIDp",&mu_plus_PIDp,"B_0_Mu_0_PIDp/D");
       NewTree->Branch("B_0_Mu_0_ProbNNe",&mu_plus_ProbNNe,"B_0_Mu_0_ProbNNe/D");
       NewTree->Branch("B_0_Mu_0_ProbNNmu",&mu_plus_ProbNNmu,"B_0_Mu_0_ProbNNmu/D");
       NewTree->Branch("B_0_Mu_0_ProbNNk",&mu_plus_ProbNNk,"B_0_Mu_0_ProbNNk/D");
       NewTree->Branch("B_0_Mu_0_ProbNNp",&mu_plus_ProbNNp,"B_0_Mu_0_ProbNNp/D");
       NewTree->Branch("B_0_Mu_0_ProbNNpi",&mu_plus_ProbNNpi,"B_0_Mu_0_ProbNNpi/D");

       NewTree->Branch("B_0_Mu_1_P",&mu_minus_P,"B_0_Mu_1_P/D");  //B_0_Mu_1
       NewTree->Branch("B_0_Mu_1_PT",&mu_minus_PT,"B_0_Mu_1_PT/D");
       NewTree->Branch("B_0_Mu_1_TRUEID",&mu_minus_TRUEID,"B_0_Mu_1_TRUEID/I");
       NewTree->Branch("B_0_Mu_1_PIDe",&mu_minus_PIDe,"B_0_Mu_1_PIDe/D");
       NewTree->Branch("B_0_Mu_1_PIDmu",&mu_minus_PIDmu,"B_0_Mu_1_PIDmu/D");
       NewTree->Branch("B_0_Mu_1_PIDK",&mu_minus_PIDK,"B_0_Mu_1_PIDK/D");
       NewTree->Branch("B_0_Mu_1_PIDp",&mu_minus_PIDp,"B_0_Mu_1_PIDp/D");
       NewTree->Branch("B_0_Mu_1_ProbNNe",&mu_minus_ProbNNe,"B_0_Mu_1_ProbNNe/D");
       NewTree->Branch("B_0_Mu_1_ProbNNmu",&mu_minus_ProbNNmu,"B_0_Mu_1_ProbNNmu/D");
       NewTree->Branch("B_0_Mu_1_ProbNNk",&mu_minus_ProbNNk,"B_0_Mu_1_ProbNNk/D");
       NewTree->Branch("B_0_Mu_1_ProbNNp",&mu_minus_ProbNNp,"B_0_Mu_1_ProbNNp/D");
       NewTree->Branch("B_0_Mu_1_ProbNNpi",&mu_minus_ProbNNpi,"B_0_Mu_1_ProbNNpi/D");
   }

   for (int evt = 0; evt < tree->GetEntries(); evt++){
       tree->GetEntry(evt);
       NewTree->Fill();
   }

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


   string dict_path = Kst2Kspiplus ? "/home/he/dgerick/Public/PID/ForDict.txt" : "/home/lhcb/kopecna/public/PID/ForDict.txt";
   ofstream myfile;
   myfile.open (dict_path,ios::app);
   myfile <<  eventTypeStr(ReferenceChannel,PHSP) << ":[\"" << year <<  "\",\"" << magnet << "\", " << NewTree->GetEntries() <<"]," << endl;

   output->Close();
   cout << "[INFO]\t\t New tree created and saved to file." << endl;
   return 1;

}
int copyForPIDAllSig(){

    for (auto& year : yearsMC(false, false,12)){
        if (copyForPID(year,"down",false,false)==0) return 0;
        if (copyForPID(year,"up",  false,false)==0) return 0;
    }

    return 1;
}
int copyForPIDAllRef(){
    for (auto& year : yearsMC(true, false,12)){
        if (copyForPID(year,"down",true,false)==0) return 0;
        if (copyForPID(year,"up",  true,false)==0) return 0;
    }
    return 1;

}

int copySameRunNumber(int year1 = 2018, int year2 = 2020, std::string mag = "down") {


    //put both magnet polarities into one TChain
    string tupleName = Kst2Kpluspi0Resolved ? "b2KstKpi0mumuResolvedTuple/DecayTree" : "b2KstKs0pimumu_Tuple/DecayTree";


    TChain * treeData1 = new TChain(tupleName.c_str());
    TChain * treeData2 = new TChain(tupleName.c_str());


    if (year1 < 2019)
        treeData1->Add(Form("%s/data/%i%s/%i*_B2Kstmumu*.root", path_to_output_KshortPiplus.c_str(), year1, mag.c_str(), year1));
    else
        treeData1->Add(Form("%s/data/%i%s/%i*_B2Kstmumu*_validation.root", Kst2Kpluspi0Resolved ? path_to_output_KplusPizero.c_str() : path_to_output_KshortPiplus.c_str(), year1, mag.c_str(), year1));

    if (year2 < 2019)
        treeData2->Add(Form("%s/data/%i%s/%i*_B2Kstmumu*.root", path_to_output_KshortPiplus.c_str(), year2, mag.c_str(), year2));
    else
        treeData2->Add(Form("%s/data/%i%s/%i*_B2Kstmumu*_validation.root", Kst2Kpluspi0Resolved ? path_to_output_KplusPizero.c_str() : path_to_output_KshortPiplus.c_str(), year2, mag.c_str(), year2));

    if(treeData1->GetEntries() == 0){
      std::cout << "[ERROR]\tNo entries found in ntuple for year " << year1 << std::endl;
      return 0;
    }
    if(treeData2->GetEntries() == 0){
      std::cout << "[ERROR]\tNo entries found in ntuple for year " << year2 << std::endl;
      return 0;
    }

    UInt_t runNumber1 = 0;
    UInt_t runNumber2 = 0;

    treeData1->SetBranchStatus("*",1);
    treeData1->SetBranchStatus("runNumber",1);
    treeData1->SetBranchAddress("runNumber",&runNumber1);

    treeData2->SetBranchStatus("*",1);
    treeData2->SetBranchStatus("runNumber",1);
    treeData2->SetBranchAddress("runNumber",&runNumber2);


    std::vector<int>vRunNumber;
    std::vector<int>vRunNumber1;
    std::vector<int>vRunNumber2;

    /*
    for (int evt1 = 0; evt1 < treeData1->GetEntries(); evt1++){
        treeData1->GetEntry(evt1);
        vRunNumber1.push_back(runNumber1);
    }
    */
    for (int evt2 = 0; evt2< treeData2->GetEntries(); evt2++){
        treeData2->GetEntry(evt2);
        vRunNumber.push_back(runNumber2);
    }

    /*
    std::sort(vRunNumber1.begin(), vRunNumber1.end());
    auto last1 = std::unique(vRunNumber1.begin(), vRunNumber1.end());
    vRunNumber1.erase(last1, vRunNumber1.end());
    for (std::vector<int>::const_iterator i = vRunNumber1.begin(); i != vRunNumber1.end(); ++i){
        cout << *i << endl;
    }
    cout << "First vector" << endl;

    std::sort(vRunNumber2.begin(), vRunNumber2.end());
    auto last2 = std::unique(vRunNumber2.begin(), vRunNumber2.end());
    vRunNumber2.erase(last2, vRunNumber2.end());

    for (std::vector<int>::const_iterator ii = vRunNumber2.begin(); ii != vRunNumber2.end(); ++ii){
        cout << *ii << endl;
    }
    cout << "Second vector" << endl;
    for (int i = 0; i < vRunNumber1.size(); i++){
        for (int ii = 0; ii < vRunNumber2.size(); ii++){
            if ( vRunNumber1[i]  == vRunNumber2[ii]){
                vRunNumber.push_back(vRunNumber1[i]);
                cout << vRunNumber1[i]  << endl;
                continue;
            }
        }
    }
    */

    std::sort(vRunNumber.begin(), vRunNumber.end());
    auto last = std::unique(vRunNumber.begin(), vRunNumber.end());
    vRunNumber.erase(last, vRunNumber.end());


    cout << "Loop to get the run numbers done" << endl;


 /*
    TTree * NewTree1 = new TTree(tupleName.c_str(),tupleName.c_str());
    for (int evt = 0; evt < treeData1->GetEntries(); evt++){
        treeData1->GetEntry(evt);
        for(unsigned rn = 0; rn < vRunNumber.size(); rn++){
            if (runNumber1 == vRunNumber.at(rn)){
                NewTree1->Fill();
                continue;
            }
        }
    }
    cout << "Loop to get the first tree done" << endl;

    TTree * NewTree2 = new TTree(tupleName.c_str(),tupleName.c_str());
    for (int evt = 0; evt < treeData2->GetEntries(); evt++){
        treeData2->GetEntry(evt);
        for(unsigned rn = 0; rn < vRunNumber.size(); rn++){
            if (runNumber2 == vRunNumber.at(rn)){
                NewTree2->Fill();
                continue;
            }
        }
        cout << evt << endl;
    }
    */

    if (vRunNumber.size() ==0) return 0;

    std::string AllCuts = "runNumber== " + to_string(vRunNumber.at(0));

    cout << "[INFO]\t\t Appling cuts:" << AllCuts << endl;

    for (vector<int>::iterator t = vRunNumber.begin()+1;
                                       t != vRunNumber.end();
                                       ++t){
        AllCuts = AllCuts + " || runNumber== " + to_string(*t) ;
        cout << "\t\t\t" << " || runNumber== "  << *t << endl;
    }

    cout << "Start copying the trees" << endl;

    TTree* NewTree1 = treeData1->CopyTree( AllCuts.c_str() );
    TTree* NewTree2 = treeData2->CopyTree( AllCuts.c_str() );

    cout << "Copying the trees done" << endl;

    TFile *output1, *output2;
    if (year1 < 2019)
        output1 = new TFile(Form("%s/data/%i%s/%i%i*_B2Kstmumu*_runNumberSame.root", path_to_output_KshortPiplus.c_str(), year1, mag.c_str(), year1, year2),"RECREATE");
    else
        output1 = new TFile(Form("%s/data/%i%s/%i%i*_B2Kstmumu*_validation_runNumberSame.root", Kst2Kpluspi0Resolved ? path_to_output_KplusPizero.c_str() : path_to_output_KshortPiplus.c_str(), year1, mag.c_str(), year1, year2),"RECREATE");


    if (year2 < 2019)
        output2 = new TFile(Form("%s/data/%i%s/%i%i*_B2Kstmumu*_runNumberSame.root", path_to_output_KshortPiplus.c_str(), year2, mag.c_str(), year1, year2),"RECREATE");
    else
        output2 = new TFile(Form("%s/data/%i%s/%i%i*_B2Kstmumu*_validation_runNumberSame.root", Kst2Kpluspi0Resolved ? path_to_output_KplusPizero.c_str() : path_to_output_KshortPiplus.c_str(), year2, mag.c_str(), year1, year2),"RECREATE");


   output1->cd();
   NewTree1->Write("",TObject::kWriteDelete);

   output2->cd();
   NewTree2->Write("",TObject::kWriteDelete);


   output1->Close();
   output2->Close();
   cout << "[INFO]\t\t New tree created and saved to file." << endl;
   return 1;

}

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

    //Load all files
    TChain* tree = 0;

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

    if(ReferenceChannel && PHSP){
        std::cout << "[WARNING]\tCannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!" << std::endl;
        PHSP = false;
    }


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

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

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


    int N = tree->GetEntries();
    coutDebug("Old files contain " + to_string(N) + " events.");
    if(N == 0){
        coutERROR("Tree has no entries. Exit!");
        return 0;
    }

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

    //activate needed branches
    tree->SetBranchStatus("mu_plus_PT",1) ;
    tree->SetBranchStatus("mu_minus_PT",1) ;
    tree->SetBranchStatus("mu_plus_ProbNNmu",1) ;
    tree->SetBranchStatus("mu_minus_ProbNNmu",1) ;
    tree->SetBranchStatus("mu_plus_PIDmu",1) ;
    tree->SetBranchStatus("mu_minus_PIDmu",1) ;
    tree->SetBranchStatus("nSPDHits",1) ;
    if (MC){
        tree->SetBranchStatus("mu_plus_TRUEPT",1) ;
        tree->SetBranchStatus("mu_minus_TRUEPT",1) ;
    }
    tree->SetBranchStatus("B_plus_L0MuonDecision_TOS",1) ;
    tree->SetBranchStatus("B_plus_L0Global_TIS",1);
    if (year > 2012) tree->SetBranchStatus("B_plus_L0DiMuonDecision_TOS",1) ;

    //open output file
    string output_path = GetInputFile(year, magnet,true,MC,ReferenceChannel,PHSP,false);
    replace(output_path,".root","_L0Muon.root");
    TFile *output = new TFile(output_path.c_str(),"RECREATE");
    output->cd();

    TTree* newTree = tree->CopyTree("B_plus_L0Global_TIS==1");

    if (newTree->GetEntries()==0) coutERROR("Failed while coping!");

    newTree->Write("",TObject::kOverwrite);

    delete newTree;
    output->Close();

    return 1;
    //TODO: delete pointers

}
int copyForL0Efficiency(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false){
    return copyForL0Efficiency(stoi(year),magnet,MC,ReferenceChannel,PHSP);
}
int copyForL0EfficiencyData(int Run){
    for (auto& year : yearsData(Run)){
        if (copyForL0Efficiency(year,"up",false, false, false)==0) return 0;
        if (copyForL0Efficiency(year,"down",false, false, false)==0) return 0;
    }
    return 1;
}
int copyForL0EfficiencyMC(int Run){
    for (auto& year : yearsMC(false,false,Run)){
        if (copyForL0Efficiency(year,"up",true, false, false)==0) return 0;
        if (copyForL0Efficiency(year,"down",true, false, false)==0) return 0;
    }
    return 1;
}
int copyForL0EfficiencyRefMC(int Run){
    for (auto& year : yearsMC(true,false,Run)){
        if (copyForL0Efficiency(year,"up",false, true, false)==0) return 0;
        if (copyForL0Efficiency(year,"down",false, true, false)==0) return 0;
    }
    return 1;
}
int copyForL0EfficiencyPHSP(int Run){
    for (auto& year : yearsMC(false,true,Run)){
        if (copyForL0Efficiency(year,"up",false, false, true)==0) return 0;
        if (copyForL0Efficiency(year,"down",false, false, true)==0) return 0;
    }
    return 1;
}
int copyForL0EfficiencyAll(int Run){
    if (copyForL0EfficiencyData(Run)==0) return 0;
    if (copyForL0EfficiencyMC(Run)==0) return 0;
    if (copyForL0EfficiencyRefMC(Run)==0) return 0;
    if (copyForL0EfficiencyPHSP(Run)==0) return 0;
    return 1;
}