
#include "../GlobalFunctions.hh"

//////////////////////////////////////////////////////////////////////
/// DTFcomparison()
/// #todo
///


string figureFormat = "eps";


std::array<string, 3> ParticlesPiZero {"B",
                                       //"muplus",
                                       //"muminus",
                                       "Kstar",
                                       "pi0",
                                       //"gamma1",
                                       //"gamma2"
                                      };

std::array<string, 3> ParticlesKShort {"B",
                                       //"muplus",
                                       //"muminus",
                                       "Kstar",
                                       "Kshort"
                                      };

std::array<string, 3> Variables{ "M",
                                 "PT",
                                 "ETA"
                               };

std::array<string,1> ExtraVariables{// "ThetaK",
    // "ThetaL",
    // "Phi",
    "OpenAngle",
};

std::array<string, 2> DTFoptions{//"DTF_PDG",
    //"PDG",
    "DTF",
    ""
};


string particleLatex(string particle){ //one doesn't have to think much when making labels and titles
    if (particle == "B")        return ("B^{+}");
    if (particle == "Bplus")    return ("B^{+}");
    if (particle == "B_plus")   return ("B^{+}");
    if (particle == "Kstar")    return ("K^{*+}");
    if (particle == "K_star")   return ("K^{*+}");
    if (particle == "Kplus")    return ("K^{+}");
    if (particle == "K_plus")   return ("K^{+}");
    if (particle == "Kshort")   return ("K_{s}");
    if (particle == "K_short")  return ("K_{s}");
    if (particle == "piplus")   return ("#pi^{+}");
    if (particle == "pi_plus")  return ("#pi^{+}");
    if (particle == "pizero")   return ("#pi^{0}");
    if (particle == "pi_zero")  return ("#pi^{0}");
    if (particle == "pi0")      return ("#pi^{0}");
    if (particle == "mu")       return ("#mu");
    if (particle == "mu_plus")  return ("#mu^{+}");
    if (particle == "muplus")   return ("#mu^{+}");
    if (particle == "mu_minus") return ("#mu^{-}");
    if (particle == "muminus")  return ("#mu^{-}");
    if (particle == "gamma")    return ("#gamma");
    cout << "Particle "<< particle <<"  not in the list!" << endl;
    return "";
}

string variableLatex(string var){
    if (var=="PT")  return ("p_{T}");
    if (var=="PX")  return ("p_{X}");
    if (var=="PY")  return ("p_{Y}");
    if (var=="PZ")  return ("p_{Z}");
    if (var=="ETA") return ("#eta");
    if (var=="M")   return ("m");
    if (var=="OpenAngle")   return ("Open angle (#gamma#gamma angle)");
    cout << "Variable "<< var <<" is not in the list!" << endl;
    return "";
}

string optionLatex(string opt){
    if (opt=="DTF_PDG") return ("#splitline{info from DTF}{with masses fixed to PDG}");
    if (opt=="DTF")     return ("info from DTF");
    if (opt=="PDG")     return ("#splitline{info from TupleTool}{with masses fixed to PDG}");
    if (opt=="")        return ("info from TupleTool");
    cout << "Option " << opt << " not in the list!" << endl;
    return "";
}

string histName(string particle, string var, string opt){
    return (variableLatex(var)+"("+ particleLatex(particle)+ "):" + optionLatex(opt));
}

void designCanvas(TCanvas* canvas){
    canvas->SetRightMargin(0.04);
    canvas->SetTopMargin(0.04);
    canvas->SetLeftMargin(0.12);
    canvas->SetBottomMargin(0.12);
}

void designPlots(TH1D *hist, int color){
    hist->SetLineWidth(2);
    hist->SetLineColor(color);
    hist->GetYaxis()->SetTitleOffset(1.7);
    hist->GetXaxis()->SetTitleOffset(1.0);
    hist->SetTitle("");
}

void designLegend(TLegend *leg, string legTitle, TH1D *hist, string opt){
    leg->SetTextSize(0.03);
    leg->SetHeader(legTitle.c_str());
    leg->AddEntry(hist,optionLatex(opt).c_str(),"l");
}

int DTFcomparison(std::string year = "2011", std::string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false, bool preselected = true, bool truthMatched = false) {
    //no point in keepeng the MC option, but it doesn't hurt neither

    TFile* output = 0;
    TChain *tree;
    if (!(MC || ReferenceChannel || PHSP)) truthMatched = false;
    if (!preselected){
        //Kst2Kpluspi0Resolved case
        if(Kst2Kpluspi0Resolved){
            if(!MC){
                tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
                if(smallSample) tree->Add(Form("%s/data/%s%s/*B2Kstmumu*13.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                else{
                    tree->Add(Form("%s/data/%s%s/*.root",path_to_data.c_str(),year.c_str(),magnet.c_str()));
                    std::cout << "Adding " << Form("%s/data/%s%s/*B2Kstmumu*.root",path_to_data.c_str(),year.c_str(),magnet.c_str()) << endl;
                }
            }
            else{
                tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
                if(ReferenceChannel){
                    if(smallSample) tree->Add(Form("%s/data/MC/RefKplusPi0/%s%s/*B2KstJpsi*1.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                    else tree->Add(Form("%s/data/MC/RefKplusPi0/%s%s/*B2KstJpsi*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                }
                else if(PHSP){
                    if(smallSample) tree->Add(Form("%s/data/PHSP/KplusPi0/%s%s/*B2Kstmumu*1.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                    else tree->Add(Form("%s/data/PHSP/KplusPi0/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                }
                else{
                    if(smallSample) tree->Add(Form("%s/data/MC/KplusPi0/%s%s/*B2Kstmumu*1.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                    else tree->Add(Form("%s/data/MC/KplusPi0/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                }
            }
        }

        //Kst2Kpluspi0Merged case
        if(Kst2Kpluspi0Merged){
            tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
            if(!MC){
                if(smallSample) tree->Add(Form("%s/data/%s%s/*B2Kstmumu*13.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                else tree->Add(Form("%s/data/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
            }
            else{
                if(smallSample) tree->Add(Form("%s/data/MC/KplusPi0/%s%s/*B2Kstmumu*1.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                else tree->Add(Form("%s/data/MC/KplusPi0/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
            }
        }

        //Kst2Kspiplus case
        if(Kst2Kspiplus){
            if(!MC){
                tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
                if(smallSample) tree->Add(Form("%s/data/%s%s/*B2Kstmumu*23.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                else tree->Add(Form("%s/data/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
            }
            else{
                tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
                if(ReferenceChannel){
                    if(smallSample) tree->Add(Form("%s/data/MC/RefKshortPiplus/%s%s/*B2KstJpsi*13.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                    else tree->Add(Form("%s/data/MC/RefKshortPiplus/%s%s/*B2KstJpsi*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                }
                else if(PHSP){
                    if(smallSample) tree->Add(Form("%s/data/PHSP/KshortPiplus/%s%s/*B2Kstmumu*1.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                    else tree->Add(Form("%s/data/PHSP/KshortPiplus/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                }
                else{
                    if(smallSample) tree->Add(Form("%s/data/MC/KshortPiplus/%s%s/*B2Kstmumu*1.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                    else tree->Add(Form("%s/data/MC/KshortPiplus/%s%s/*B2Kstmumu*.root",path_to_data.c_str(), year.c_str(),magnet.c_str()));
                }
            }
        }

    }
    else{ //preselected
        string name = (truthMatched ? "DecayTreeTruthMatched" : "DecayTree");
        tree=new TChain(name.c_str());
        // Kst2Kpluspi0Resolved case
        if(Kst2Kpluspi0Resolved){
            if(!MC){
                tree->Add(Form("%s/data/%s%s/%s%s_pi0Resolved.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                cout << Form("%s/data/%s%s/%s%s_pi0Resolved.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()) << endl;
            }
            else{
                if(ReferenceChannel){
                    tree->Add(Form("%s/data/MC/RefKplusPi0/%s%s/%s%s_pi0Resolved.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
                else if(PHSP){
                    tree->Add(Form("%s/data/PHSP/KplusPi0/%s%s/%s%s_pi0Resolved.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
                else{
                    tree->Add(Form("%s/data/MC/KplusPi0/%s%s/%s%s_pi0Resolved.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
            }
        }
        // Kst2Kpluspi0Merged case
        if(Kst2Kpluspi0Merged){
            if(!MC){
                tree->Add(Form("%s/data/%s%s/%s%s_pi0Merged.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
            }
            else{
                if(ReferenceChannel){
                    tree->Add(Form("%s/data/MC/RefKplusPi0/%s%s/%s%s_pi0Merged.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
                else if(PHSP){
                    tree->Add(Form("%s/data/PHSP/KplusPi0/%s%s/%s%s_pi0Merged.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
                else{
                    tree->Add(Form("%s/data/MC/KplusPi0/%s%s/%s%s_pi0Mergedroot",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
            }
        }
        // Kst2Kspiplus case
        if(Kst2Kspiplus){
            if(!MC){
                tree->Add(Form("%s/data/%s%s/%s%s_piplus.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
            }
            else{
                if(ReferenceChannel){
                    tree->Add(Form("%s/data/MC/RefKshortPiplus/%s%s/%s%s_piplus.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
                else if(PHSP){
                    tree->Add(Form("%s/data/PHSP/KshortPiplus/%s%s/%s%s_piplus.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
                else{
                    tree->Add(Form("%s/data/MC/KshortPiplus/%s%s/%s%s_piplus.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()));
                }
            }
        }
    }

    int N = tree->GetEntries();

    if(N == 0){
	std::cout << "[ERROR]\t\tNo events found!" << std::endl;
	return 0;
    }

    // MASSES HISTOGRAMS
    string tmpName = "";


    //pi0
    tmpName = "m_{"+ particleLatex("pi0") +"}[MeV]";
    TH1 * h_pi0_M_DTF_PDG= new TH1D("h_pi0_M_DTF_PDG", histName("pi0","M", "DTF_PDG").c_str(), 70, 100, 170);
    h_pi0_M_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_M_DTF_PDG->GetYaxis()->SetTitle("Counts / 1MeV");
    TH1 * h_pi0_M_DTF= new TH1D("h_pi0_M_DTF" , histName("pi0","M", "DTF").c_str(), 70, 100, 170);
    h_pi0_M_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_M_DTF->GetYaxis()->SetTitle("Counts / 1MeV");

    tmpName = "m_{"+ particleLatex("pizero") +"}[MeV]";
    TH1 * h_pi0_M_PDG= new TH1D("h_pi0_M_PDG" , histName("pi0","M", "PDG").c_str(), 70, 100, 170);
    h_pi0_M_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_M_PDG->GetYaxis()->SetTitle("Counts / 1MeV");
    TH1 * h_pi0_M= new TH1D("h_pi0_M" , histName("pi0","M", "").c_str(), 70, 100, 170);
    h_pi0_M->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_M->GetYaxis()->SetTitle("Counts / 1MeV");

    //Ks
    tmpName = "m_{" + particleLatex("Kshort") + "}[MeV]";
    TH1 * h_Kshort_M_DTF_PDG= new TH1D("h_Kshort_M_DTF_PDG" , histName("Kshort","M", "DTF_PDG").c_str(), 100, 400, 600);
    h_Kshort_M_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_M_DTF_PDG->GetYaxis()->SetTitle("Counts / 2MeV");
    TH1 * h_Kshort_M_DTF= new TH1D("h_Kshort_M_DTF" , histName("Kshort","M", "DTF").c_str(), 100, 400, 600);
    h_Kshort_M_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_M_DTF->GetYaxis()->SetTitle("Counts / 2MeV");
    TH1 * h_Kshort_M_PDG= new TH1D("h_Kshort_M_PDG" ,histName("Kshort","M", "PDG").c_str(), 100, 400, 600);
    h_Kshort_M_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_M_PDG->GetYaxis()->SetTitle("Counts / 2MeV");
    TH1 * h_Kshort_M= new TH1D("h_Kshort_M" , histName("Kshort","M", "").c_str(), 100, 400, 600);
    h_Kshort_M->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_M->GetYaxis()->SetTitle("Counts / 2MeV");

    //K*
    if (Kst2Kpluspi0Resolved) tmpName = "m_{" + particleLatex("Kplus") + particleLatex("pi0") + "}[MeV]";
    else if (Kst2Kspiplus) tmpName = "m_{" + particleLatex("Kshort") + particleLatex("piplus") + "}[MeV]";
    TH1 * h_Kstar_M_DTF_PDG= new TH1D("h_Kstar_M_DTF_PDG" ,histName("Kstar","M", "DTF_PDG").c_str(), 100, 700, 1100);
    h_Kstar_M_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_M_DTF_PDG->GetYaxis()->SetTitle("Counts / 4MeV");
    TH1 * h_Kstar_M_DTF= new TH1D("h_Kstar_M_DTF" , histName("Kstar","M", "DTF").c_str(), 100, 700, 1100);
    h_Kstar_M_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_M_DTF->GetYaxis()->SetTitle("Counts / 4MeV");
    if (Kst2Kpluspi0Resolved) tmpName = "m_{" + particleLatex("Kplus") + particleLatex("pizero") + "}[MeV]";
    TH1 * h_Kstar_M_PDG= new TH1D("h_Kstar_M_PDG" , histName("Kstar","M", "PDG").c_str(), 100, 700, 1100);
    h_Kstar_M_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_M_PDG->GetYaxis()->SetTitle("Counts / 4MeV");
    TH1 * h_Kstar_M= new TH1D("h_Kstar_M" , histName("Kstar","M", "").c_str(), 100, 700, 1100);
    h_Kstar_M->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_M->GetYaxis()->SetTitle("Counts / 4MeV");

    //B+
    if (Kst2Kpluspi0Resolved) tmpName = "m_{" + particleLatex("Kplus") + particleLatex("pi0") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    else if (Kst2Kspiplus) tmpName = "m_{" + particleLatex("Kshort") + particleLatex("piplus") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    TH1 * h_B_M_DTF_PDG= new TH1D("h_B_M_DTF_PDG" , histName("Bplus","M", "DTF_PDG").c_str(), (7000-cut_B_plus_M_low)/20, cut_B_plus_M_low, 7000);
    h_B_M_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_M_DTF_PDG->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_B_M_DTF= new TH1D("h_B_M_DTF" , histName("Bplus","M", "DTF").c_str(), (7000-cut_B_plus_M_low)/20, cut_B_plus_M_low, 7000);
    h_B_M_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_M_DTF->GetYaxis()->SetTitle("Counts / 20MeV");

    if (Kst2Kpluspi0Resolved) tmpName = "m_{" + particleLatex("Kplus") + particleLatex("pizero") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    TH1 * h_B_M_PDG= new TH1D("h_B_M_PDG", histName("Bplus","M", "PDG").c_str(), (7000-cut_B_plus_M_low)/20, cut_B_plus_M_low, 7000);
    h_B_M_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_M_PDG->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_B_M= new TH1D("h_B_M" , histName("Bplus","M", "").c_str(), (7000-cut_B_plus_M_low)/20, cut_B_plus_M_low, 7000);
    h_B_M->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_M->GetYaxis()->SetTitle("Counts / 20MeV");


    //pT HISTOGRAMS
    //pi0
    tmpName = variableLatex("PT") + "{"+ particleLatex("pi0") +"}[MeV]";
    TH1 * h_pi0_PT_DTF_PDG= new TH1D("h_pi0_PT_DTF_PDG" , histName("pi0","PT", "DTF_PDG").c_str(),  100, 0, 5000);
    h_pi0_PT_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_PT_DTF_PDG->GetYaxis()->SetTitle("Counts / 50MeV");
    TH1 * h_pi0_PT_DTF= new TH1D("h_pi0_PT_DTF" ,  histName("pizero","PT", "DTF").c_str(),100, 0, 5000);
    h_pi0_PT_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_PT_DTF->GetYaxis()->SetTitle("Counts / 50MeV");
    TH1 * h_pi0_PT_PDG= new TH1D("h_pi0_PT_PDG" ,  histName("pi0","PT", "PDG").c_str(), 100, 0, 5000);;
    h_pi0_PT_PDG->GetXaxis()->SetTitle("p_{T}{#pi^{0}[MeV]");
    h_pi0_PT_PDG->GetYaxis()->SetTitle("Counts / 50MeV");
    TH1 * h_pi0_PT= new TH1D("h_pi0_PT" , histName("pi0","PT", "").c_str(),100, 0, 5000);
    h_pi0_PT->GetXaxis()->SetTitle("p_{T}{#pi^{0}[MeV]");
    h_pi0_PT->GetYaxis()->SetTitle("Counts / 50MeV");

    //Ks
    tmpName = variableLatex("PT") + "{" + particleLatex("Kshort") + "}[MeV]";
    TH1 * h_Kshort_PT_DTF_PDG= new TH1D("h_Kshort_PT_DTF_PDG" , histName("Kshort","PT", "DTF_PDG").c_str(), 100, 0, 15000);
    h_Kshort_PT_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_PT_DTF_PDG->GetYaxis()->SetTitle("Counts / 150MeV");
    TH1 * h_Kshort_PT_DTF= new TH1D("h_Kshort_PT_DTF" ,  histName("Kshort","PT", "DTF").c_str(), 100, 0, 15000);
    h_Kshort_PT_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_PT_DTF->GetYaxis()->SetTitle("Counts / 150MeV");
    TH1 * h_Kshort_PT_PDG= new TH1D("h_Kshort_PT_PDG" ,  histName("Kshort","PT", "PDG").c_str(), 100,  0, 15000);
    h_Kshort_PT_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_PT_PDG->GetYaxis()->SetTitle("Counts / 150MeV");
    TH1 * h_Kshort_PT= new TH1D("h_Kshort_PT" ,  histName("Kshort","PT", "").c_str(), 100, 0, 15000);
    h_Kshort_PT->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_PT->GetYaxis()->SetTitle("Counts / 150MeV");

    //K*
    if (Kst2Kpluspi0Resolved) tmpName = variableLatex("PT") + "{" + particleLatex("Kplus") + particleLatex("pi0") + "}[MeV]";
    else if (Kst2Kspiplus) tmpName = variableLatex("PT") + "{" + particleLatex("Kshort") + particleLatex("piplus") + "}[MeV]";
    TH1 * h_Kstar_PT_DTF_PDG= new TH1D("h_Kstar_PT_DTF_PDG" ,  histName("Kstar","PT", "DTF_PDG").c_str(), 100, 0, 15000);
    h_Kstar_PT_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_PT_DTF_PDG->GetYaxis()->SetTitle("Counts / 150MeV");
    TH1 * h_Kstar_PT_DTF= new TH1D("h_Kstar_PT_DTF" ,  histName("Kstar","PT", "DTF").c_str(), 100, 0, 15000);
    h_Kstar_PT_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_PT_DTF->GetYaxis()->SetTitle("Counts / 150MeV");
    TH1 * h_Kstar_PT_PDG= new TH1D("h_Kstar_PT_PDG" , histName("Kstar","PT", "PDG").c_str(),100, 0, 15000);
    h_Kstar_PT_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_PT_PDG->GetYaxis()->SetTitle("Counts / 150MeV");
    TH1 * h_Kstar_PT= new TH1D("h_Kstar_PT" ,  histName("Kstar","PT", "").c_str(), 100, 0, 15000);
    h_Kstar_PT->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_PT->GetYaxis()->SetTitle("Counts / 150MeV");

    //B+
    if (Kst2Kpluspi0Resolved) tmpName = variableLatex("PT") + "{" + particleLatex("Kplus") + particleLatex("pi0") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    else if (Kst2Kspiplus) tmpName = variableLatex("PT") + "{" + particleLatex("Kshort") + particleLatex("piplus") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    TH1 * h_B_PT_DTF_PDG= new TH1D("h_B_PT_DTF_PDG" , histName("Bplus","PT", "DTF_PDG").c_str(), 100, 0, 40000);
    h_B_PT_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_PT_DTF_PDG->GetYaxis()->SetTitle("Counts / 400MeV");
    TH1 * h_B_PT_DTF= new TH1D("h_B_PT_DTF" , histName("Bplus","PT", "DTF").c_str(), 100, 0, 40000);
    h_B_PT_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_PT_DTF->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_B_PT_PDG= new TH1D("h_B_PT_PDG" , histName("Bplus","PT", "PDG").c_str(), 100, 0, 40000);
    h_B_PT_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_PT_PDG->GetYaxis()->SetTitle("Counts / 20MeV");
    TH1 * h_B_PT= new TH1D("h_B_PT" ,  histName("Bplus","PT", "").c_str(), 100,  0, 40000);
    h_B_PT->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_PT->GetYaxis()->SetTitle("Counts / 20MeV");

    //eta HISTOGRAMS #TODO: hist titles
    //pi0
    tmpName = variableLatex("ETA") + "{"+ particleLatex("pi0") +"}[MeV]";
    TH1 * h_pi0_ETA_DTF_PDG= new TH1D("h_pi0_ETA_DTF_PDG" , histName("pi0","ETA", "DTF_PDG").c_str(), 35, 1.5, 5);
    h_pi0_ETA_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_ETA_DTF_PDG->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_pi0_ETA_DTF= new TH1D("h_pi0_ETA_DTF" , histName("pi0","ETA", "DTF").c_str(),  35, 1.5, 5);
    h_pi0_ETA_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_ETA_DTF->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_pi0_ETA_PDG= new TH1D("h_pi0_ETA_PDG" , histName("pi0","ETA", "PDG").c_str(),  35, 1.5, 5);
    h_pi0_ETA_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_ETA_PDG->GetYaxis()->SetTitle("Counts /  0.1");
    TH1 * h_pi0_ETA= new TH1D("h_pi0_ETA" , histName("pi0","ETA", "").c_str(),  35, 1.5, 5);
    h_pi0_ETA->GetXaxis()->SetTitle(tmpName.c_str());
    h_pi0_ETA->GetYaxis()->SetTitle("Counts / 0.1");

    //Ks
    tmpName = variableLatex("ETA") + "{" + particleLatex("Kshort") + "}[MeV]";
    TH1 * h_Kshort_ETA_DTF_PDG= new TH1D("h_Kshort_ETA_DTF_PDG" , histName("Kshort","ETA", "DTF_PDG").c_str(), 35, 1.5, 5);
    h_Kshort_ETA_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_ETA_DTF_PDG->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_Kshort_ETA_DTF= new TH1D("h_Kshort_ETA_DTF" , histName("Kshort","ETA", "DTF").c_str(), 35, 1.5, 5);
    h_Kshort_ETA_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_ETA_DTF->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_Kshort_ETA_PDG= new TH1D("h_Kshort_ETA_PDG" , histName("Kshort","ETA", "PDG").c_str(), 35, 1.5, 5);
    h_Kshort_ETA_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_ETA_PDG->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_Kshort_ETA= new TH1D("h_Kshort_ETA" ,  histName("Kshort","ETA", "").c_str(), 35, 1.5, 5);
    h_Kshort_ETA->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kshort_ETA->GetYaxis()->SetTitle("Counts / 0.1");

    //K*
    if (Kst2Kpluspi0Resolved) tmpName = variableLatex("ETA") + "{" + particleLatex("Kplus") + particleLatex("pi0") + "}[MeV]";
    else if (Kst2Kspiplus) tmpName = variableLatex("ETA") + "{" + particleLatex("Kshort") + particleLatex("piplus") + "}[MeV]";
    TH1 * h_Kstar_ETA_DTF_PDG= new TH1D("h_Kstar_ETA_DTF_PDG" , histName("Kstar","ETA", "DTF_PDG").c_str(), 35, 1.5, 5);
    h_Kstar_ETA_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_ETA_DTF_PDG->GetYaxis()->SetTitle("Counts /  0.1");
    TH1 * h_Kstar_ETA_DTF= new TH1D("h_Kstar_ETA_DTF" , histName("Kstar","ETA", "DTF").c_str(), 35, 1.5, 5);
    h_Kstar_ETA_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_ETA_DTF->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_Kstar_ETA_PDG= new TH1D("h_Kstar_ETA_PDG" , histName("Kstar","ETA", "PDG").c_str(), 35, 1.5, 5);
    h_Kstar_ETA_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_ETA_PDG->GetYaxis()->SetTitle("Counts /  0.1");
    TH1 * h_Kstar_ETA= new TH1D("h_Kstar_ETA" , histName("Kstar","ETA", "").c_str(), 35, 1.5, 5);
    h_Kstar_ETA->GetXaxis()->SetTitle(tmpName.c_str());
    h_Kstar_ETA->GetYaxis()->SetTitle("Counts / 0.1");

    //B+
    if (Kst2Kpluspi0Resolved) tmpName = variableLatex("ETA") + "{" + particleLatex("Kplus") + particleLatex("pi0") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    else if (Kst2Kspiplus) tmpName = variableLatex("ETA") + "{" + particleLatex("Kshort") + particleLatex("piplus") + particleLatex("mu") +  particleLatex("mu") + "}[MeV]";
    TH1 * h_B_ETA_DTF_PDG= new TH1D("h_B_ETA_DTF_PDG" , histName("Bplus","ETA", "DTF_PDG").c_str(), 45, 1.5, 6);
    h_B_ETA_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_ETA_DTF_PDG->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_B_ETA_DTF= new TH1D("h_B_ETA_DTF" , histName("Bplus","ETA", "DTF").c_str(), 45, 1.5, 6);
    h_B_ETA_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_ETA_DTF->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_B_ETA_PDG= new TH1D("h_B_ETA_PDG" , histName("Bplus","ETA", "PDG").c_str(),  45, 1.5, 6);
    h_B_ETA_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_ETA_PDG->GetYaxis()->SetTitle("Counts / 0.1");
    TH1 * h_B_ETA = new TH1D("h_B_ETA" , histName("Bplus","ETA", "").c_str(), 45, 1.5, 6);
    h_B_ETA->GetXaxis()->SetTitle(tmpName.c_str());
    h_B_ETA->GetYaxis()->SetTitle("Counts / 0.1");

    //Angular description
    TH1 * h_ctk = new TH1D("h_ctk" ,"#theta_{k}", 73, -3.15, 3.15);
    h_ctk->GetXaxis()->SetTitle("#theta_{k}");
    h_ctk->GetYaxis()->SetTitle("Counts /  0.1");
    TH1 * h_ctl = new TH1D("h_ctl" , "#theta_{l}", 73, -3.15, 3.15);
    h_ctl->GetXaxis()->SetTitle("#theta_{l}");
    h_ctl->GetYaxis()->SetTitle("Counts /  0.1");
    TH1 * h_sp = new TH1D("h_sp" ,"#phi", 73, -3.15, 3.15);
    h_sp->GetXaxis()->SetTitle("#phi");
    h_sp->GetYaxis()->SetTitle("Counts /  0.1");

    //Open Angle
    tmpName ="Open Angle";
    TH1 * h_OpenAngle_DTF_PDG= new TH1D("h_OpenAngle_DTF_PDG" , ("Open Angle" + optionLatex("DTF_PDG")).c_str(),  40, 0, 0.2);
    h_OpenAngle_DTF_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_OpenAngle_DTF_PDG->GetYaxis()->SetTitle("Counts / 0.005");
    TH1 * h_OpenAngle_DTF= new TH1D("h_OpenAngle_DTF" ,  ("Open Angle" + optionLatex("DTF")).c_str(), 40, 0, 0.2);
    h_OpenAngle_DTF->GetXaxis()->SetTitle(tmpName.c_str());
    h_OpenAngle_DTF->GetYaxis()->SetTitle("Counts / 0.005");
    TH1 * h_OpenAngle_PDG= new TH1D("h_OpenAngle_PDG" , ("Open Angle" + optionLatex("PDG")).c_str(), 40, 0, 0.2);
    h_OpenAngle_PDG->GetXaxis()->SetTitle(tmpName.c_str());
    h_OpenAngle_PDG->GetYaxis()->SetTitle("Counts / 0.005");
    TH1 * h_OpenAngle = new TH1D("h_OpenAngle" ,  ("Open Angle" + optionLatex("")).c_str(),  40, 0, 0.2);
    h_OpenAngle->GetXaxis()->SetTitle(tmpName.c_str());
    h_OpenAngle->GetYaxis()->SetTitle("Counts / 0.05");


    /////////////////////
    // Detector values //
    /////////////////////
    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;

    //Other variables:
    Double_t OpenAngle;

    /////////////////////
    //   DTF values    //
    /////////////////////
    TLorentzVector LorVec_DTF_B_plus;
    TLorentzVector LorVec_DTF_K_star_plus;
    TLorentzVector LorVec_DTF_K_plus;
    TLorentzVector LorVec_DTF_pi_zero_merged;
    TLorentzVector LorVec_DTF_pi_zero_resolved;
    TLorentzVector LorVec_DTF_K_short;
    TLorentzVector LorVec_DTF_pi_plus;
    TLorentzVector LorVec_DTF_KS_pi_plus;
    TLorentzVector LorVec_DTF_KS_pi_minus;
    TLorentzVector LorVec_DTF_mu_minus;
    TLorentzVector LorVec_DTF_mu_plus;
    TLorentzVector LorVec_DTF_gamma1;
    TLorentzVector LorVec_DTF_gamma2;
    TLorentzVector LorVec_DTF_DiGamma;

    //potentially misIDed particles:
    TLorentzVector LorVec_DTF_pi_plus_MisIDed;
    TLorentzVector LorVec_DTF_mu_plus_MisIDed;
    TLorentzVector LorVec_DTF_B_plus_From_DoubleMisIDed;
    TLorentzVector LorVec_DTF_Jpsi_From_DoubleMisIDed;
    TLorentzVector LorVec_DTF_K_star_plus_From_DoubleMisIDed;

    //recombined particles:
    TLorentzVector LorVec_DTF_DiMuon;
    TLorentzVector LorVec_DTF_KaonPion;
    TLorentzVector LorVec_DTF_ReconB;
    TLorentzVector LorVec_DTF_ReconBPiZero;
    TLorentzVector LorVec_DTF_KshortMuMu;
    TLorentzVector LorVec_DTF_KplusMuMu;

    //Other variables:
    Double_t OpenAngle_DTF;

    /////////////////////
    //  PDG+DTF values //
    /////////////////////
    TLorentzVector LorVec_DTF_PDG_B_plus;
    TLorentzVector LorVec_DTF_PDG_K_star_plus;
    TLorentzVector LorVec_DTF_PDG_K_plus;
    TLorentzVector LorVec_DTF_PDG_pi_zero_merged;
    TLorentzVector LorVec_DTF_PDG_pi_zero_resolved;
    TLorentzVector LorVec_DTF_PDG_K_short;
    TLorentzVector LorVec_DTF_PDG_pi_plus;
    TLorentzVector LorVec_DTF_PDG_KS_pi_plus;
    TLorentzVector LorVec_DTF_PDG_KS_pi_minus;
    TLorentzVector LorVec_DTF_PDG_mu_minus;
    TLorentzVector LorVec_DTF_PDG_mu_plus;
    TLorentzVector LorVec_DTF_PDG_gamma1;
    TLorentzVector LorVec_DTF_PDG_gamma2;
    TLorentzVector LorVec_DTF_PDG_DiGamma;

    //potentially misIDed particles:
    TLorentzVector LorVec_DTF_PDG_pi_plus_MisIDed;
    TLorentzVector LorVec_DTF_PDG_mu_plus_MisIDed;
    TLorentzVector LorVec_DTF_PDG_B_plus_From_DoubleMisIDed;
    TLorentzVector LorVec_DTF_PDG_Jpsi_From_DoubleMisIDed;
    TLorentzVector LorVec_DTF_PDG_K_star_plus_From_DoubleMisIDed;

    //recombined particles:
    TLorentzVector LorVec_DTF_PDG_DiMuon;
    TLorentzVector LorVec_DTF_PDG_KaonPion;
    TLorentzVector LorVec_DTF_PDG_ReconB;
    TLorentzVector LorVec_DTF_PDG_ReconBPiZero;
    TLorentzVector LorVec_DTF_PDG_KshortMuMu;
    TLorentzVector LorVec_DTF_PDG_KplusMuMu;

    //Other variables:
    Double_t OpenAngle_DTF_PDG;

    /////////////////////
    //   PDG values    //
    /////////////////////
    TLorentzVector LorVec_PDG_B_plus;
    TLorentzVector LorVec_PDG_K_star_plus;
    TLorentzVector LorVec_PDG_K_plus;
    TLorentzVector LorVec_PDG_pi_zero_merged;
    TLorentzVector LorVec_PDG_pi_zero_resolved;
    TLorentzVector LorVec_PDG_K_short;
    TLorentzVector LorVec_PDG_pi_plus;
    TLorentzVector LorVec_PDG_KS_pi_plus;
    TLorentzVector LorVec_PDG_KS_pi_minus;
    TLorentzVector LorVec_PDG_mu_minus;
    TLorentzVector LorVec_PDG_mu_plus;
    TLorentzVector LorVec_PDG_gamma1;
    TLorentzVector LorVec_PDG_gamma2;
    TLorentzVector LorVec_PDG_DiGamma;

    //potentially misIDed particles:
    TLorentzVector LorVec_PDG_pi_plus_MisIDed;
    TLorentzVector LorVec_PDG_mu_plus_MisIDed;
    TLorentzVector LorVec_PDG_B_plus_From_DoubleMisIDed;
    TLorentzVector LorVec_PDG_Jpsi_From_DoubleMisIDed;
    TLorentzVector LorVec_PDG_K_star_plus_From_DoubleMisIDed;

    //recombined particles:
    TLorentzVector LorVec_PDG_DiMuon;
    TLorentzVector LorVec_PDG_KaonPion;
    TLorentzVector LorVec_PDG_ReconB;
    TLorentzVector LorVec_PDG_ReconBPiZero;
    TLorentzVector LorVec_PDG_KshortMuMu;
    TLorentzVector LorVec_PDG_KplusMuMu;

    //Other variables:
    Double_t OpenAngle_PDG;

    /////////////////////
    //    Branches     //
    /////////////////////

    tree->SetBranchStatus("*",0);
    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( "*Phi",1 );
    tree->SetBranchStatus( "*ThetaL",1 );
    tree->SetBranchStatus( "*ThetaK",1 );

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


    ///K*+
    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.;

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

    Float_t K_plus_PX_DTF[100];
    Float_t K_plus_PY_DTF[100];
    Float_t K_plus_PZ_DTF[100];
    Float_t K_plus_PE_DTF[100];

    ///pi0 (merged)
    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.;

    ///pi0 (resolved)
    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.;

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

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

    Float_t Ks_pi_plus_PX_DTF[100];
    Float_t Ks_pi_plus_PY_DTF[100];
    Float_t Ks_pi_plus_PZ_DTF[100];
    Float_t Ks_pi_plus_PE_DTF[100];

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

    Float_t Ks_pi_minus_PX_DTF[100];
    Float_t Ks_pi_minus_PY_DTF[100];
    Float_t Ks_pi_minus_PZ_DTF[100];
    Float_t Ks_pi_minus_PE_DTF[100];

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

    Float_t pi_plus_PX_DTF[100];
    Float_t pi_plus_PY_DTF[100];
    Float_t pi_plus_PZ_DTF[100];
    Float_t pi_plus_PE_DTF[100];

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

    Float_t mu_minus_PX_DTF[100];
    Float_t mu_minus_PY_DTF[100];
    Float_t mu_minus_PZ_DTF[100];
    Float_t mu_minus_PE_DTF[100];

    ///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 mu_plus_PX_DTF[100];
    Float_t mu_plus_PY_DTF[100];
    Float_t mu_plus_PZ_DTF[100];
    Float_t mu_plus_PE_DTF[100];

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

    Float_t gamma1_PX_DTF[100];
    Float_t gamma1_PY_DTF[100];
    Float_t gamma1_PZ_DTF[100];
    Float_t gamma1_PE_DTF[100];

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

    Float_t gamma2_PX_DTF[100];
    Float_t gamma2_PY_DTF[100];
    Float_t gamma2_PZ_DTF[100];
    Float_t gamma2_PE_DTF[100];

    //angles
    Double_t ThetaL = 0.;
    Double_t ThetaK = 0.;
    Double_t Phi = 0.;



    //particle masses and momenta
    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( "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 );

    if(Kst2Kpluspi0Merged){ //TODO: in case one really wants to do that, add DTF info
        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 );
    }
    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( "gamma1_PX" , &gamma1_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PX" , &gamma1_PX_DTF );
        tree -> SetBranchAddress( "gamma1_PY" , &gamma1_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PY" , &gamma1_PY_DTF );
        tree -> SetBranchAddress( "gamma1_PZ" , &gamma1_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PZ" , &gamma1_PZ_DTF );
        tree -> SetBranchAddress( "gamma1_PE" , &gamma1_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PE" , &gamma1_PE_DTF );
        tree -> SetBranchAddress( "gamma2_PX" , &gamma2_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PX" , &gamma2_PX_DTF );
        tree -> SetBranchAddress( "gamma2_PY" , &gamma2_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PY" , &gamma2_PY_DTF );
        tree -> SetBranchAddress( "gamma2_PZ" , &gamma2_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PZ" , &gamma2_PZ_DTF );
        tree -> SetBranchAddress( "gamma2_PE" , &gamma2_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PE" , &gamma2_PE_DTF );

    }
    if(Kst2Kpluspi0Merged || Kst2Kpluspi0Resolved){
        tree -> SetBranchAddress( "K_plus_PX" , &K_plus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PX" , &K_plus_PX_DTF );
        tree -> SetBranchAddress( "K_plus_PY" , &K_plus_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PY" , &K_plus_PY_DTF );
        tree -> SetBranchAddress( "K_plus_PZ" , &K_plus_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PZ" , &K_plus_PZ_DTF );
        tree -> SetBranchAddress( "K_plus_PE" , &K_plus_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PE" , &K_plus_PE_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( "pi_plus_PX" , &pi_plus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PX" , &pi_plus_PX_DTF );
        tree -> SetBranchAddress( "pi_plus_PY" , &pi_plus_PY );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PY" , &pi_plus_PY_DTF );
        tree -> SetBranchAddress( "pi_plus_PZ" , &pi_plus_PZ );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PZ" , &pi_plus_PZ_DTF );
        tree -> SetBranchAddress( "pi_plus_PE" , &pi_plus_PE );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PE" , &pi_plus_PE_DTF );

        tree -> SetBranchAddress( "Ks_pi_minus_PX" , &Ks_pi_minus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PX" , &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" , &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" , &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" , &Ks_pi_minus_PE_DTF );

        tree -> SetBranchAddress( "Ks_pi_plus_PX" , &Ks_pi_plus_PX );
        tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PX" , &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" , &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" , &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" , &Ks_pi_plus_PE_DTF );
    }
    tree -> SetBranchAddress( "mu_minus_PX" , &mu_minus_PX );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PX" , &mu_minus_PX_DTF );
    tree -> SetBranchAddress( "mu_minus_PY" , &mu_minus_PY );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PY" , &mu_minus_PY_DTF );
    tree -> SetBranchAddress( "mu_minus_PZ" , &mu_minus_PZ );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PZ" , &mu_minus_PZ_DTF );
    tree -> SetBranchAddress( "mu_minus_PE" , &mu_minus_PE );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PE" , &mu_minus_PE_DTF );

    tree -> SetBranchAddress( "mu_plus_PX" , &mu_plus_PX );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PX" , &mu_plus_PX_DTF );
    tree -> SetBranchAddress( "mu_plus_PY" , &mu_plus_PY );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PY" , &mu_plus_PY_DTF );
    tree -> SetBranchAddress( "mu_plus_PZ" , &mu_plus_PZ );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PZ" , &mu_plus_PZ_DTF );
    tree -> SetBranchAddress( "mu_plus_PE" , &mu_plus_PE );
    tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PE" , &mu_plus_PE_DTF );


    tree -> SetBranchAddress( "B_plus_ThetaL" , &ThetaL);
    tree -> SetBranchAddress( "B_plus_ThetaK" , &ThetaK);
    tree -> SetBranchAddress( "B_plus_Phi" , &Phi );


    //////////////////////////////////
    /////                       //////
    ////     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 + 1 << "/" << N
                 << endl;
        }

        tree->GetEntry(i);

        /////////////////////
        //  PDG+DTF values //
        /////////////////////
        LorVec_DTF_PDG_mu_minus.SetXYZM         (mu_minus_PX_DTF[0],           mu_minus_PY_DTF[0],            mu_minus_PZ_DTF[0],            PDGMASS.MU);
        LorVec_DTF_PDG_mu_plus.SetXYZM          (mu_plus_PX_DTF[0],            mu_plus_PY_DTF[0],             mu_plus_PZ_DTF[0],             PDGMASS.MU);
        if(Kst2Kspiplus){
            /// Ks, pi+
            LorVec_DTF_PDG_pi_plus.SetXYZM          (pi_plus_PX_DTF[0],            pi_plus_PY_DTF[0],             pi_plus_PZ_DTF[0],             PDGMASS.PI_PLUS);
            LorVec_DTF_PDG_KS_pi_plus.SetXYZM       (Ks_pi_plus_PX_DTF[0],         Ks_pi_plus_PY_DTF[0],          Ks_pi_plus_PZ_DTF[0],          PDGMASS.PI_PLUS);
            LorVec_DTF_PDG_KS_pi_minus.SetXYZM      (Ks_pi_minus_PX_DTF[0],        Ks_pi_minus_PY_DTF[0],         Ks_pi_minus_PZ_DTF[0],         PDGMASS.PI_PLUS);
            LorVec_DTF_PDG_pi_plus_MisIDed.SetXYZM  (pi_plus_PX_DTF[0],            pi_plus_PY_DTF[0],             pi_plus_PZ_DTF[0],             PDGMASS.MU);
            LorVec_DTF_PDG_mu_plus_MisIDed.SetXYZM  (mu_plus_PX_DTF[0],            mu_plus_PY_DTF[0],             mu_plus_PZ_DTF[0],             PDGMASS.PI_PLUS);
        }
        else{
            /// K+, pi0
            LorVec_DTF_PDG_K_plus.SetXYZM           (K_plus_PX_DTF[0],             K_plus_PY_DTF[0],              K_plus_PZ_DTF[0],              PDGMASS.K_PLUS);
            LorVec_DTF_PDG_gamma1.SetXYZM           (gamma1_PX_DTF[0],             gamma1_PY_DTF[0],              gamma1_PZ_DTF[0],              PDGMASS.GAMMA);
            LorVec_DTF_PDG_gamma2.SetXYZM           (gamma2_PX_DTF[0],             gamma2_PY_DTF[0],              gamma2_PZ_DTF[0],              PDGMASS.GAMMA);
        }
        //directly recombine the non-stable particles for later use (K*+ not really needed/used for Kst2Kspiplus)
        if(Kst2Kspiplus){
            LorVec_DTF_PDG_K_short = LorVec_DTF_PDG_KS_pi_plus + LorVec_DTF_PDG_KS_pi_minus;
            LorVec_DTF_PDG_KaonPion = LorVec_DTF_PDG_K_short + LorVec_DTF_PDG_pi_plus;
        }
        if(Kst2Kpluspi0Resolved){
            LorVec_DTF_PDG_DiGamma = LorVec_DTF_PDG_gamma1 + LorVec_DTF_PDG_gamma2;
            LorVec_DTF_PDG_pi_zero_resolved = LorVec_DTF_PDG_DiGamma;
            LorVec_DTF_PDG_KaonPion = LorVec_DTF_PDG_K_plus + LorVec_DTF_PDG_pi_zero_resolved;
        }
        LorVec_DTF_PDG_K_star_plus = LorVec_DTF_PDG_KaonPion;
	LorVec_DTF_PDG_DiMuon = LorVec_DTF_PDG_mu_minus     + LorVec_DTF_PDG_mu_plus;
        LorVec_DTF_PDG_B_plus = LorVec_DTF_PDG_K_star_plus  + LorVec_DTF_PDG_DiMuon;

        /////////////////////
        //  DTF-values     //
        /////////////////////
        LorVec_DTF_mu_minus = TLorentzVector        (mu_minus_PX_DTF[0],           mu_minus_PY_DTF[0],            mu_minus_PZ_DTF[0],            mu_minus_PE_DTF[0]);
        LorVec_DTF_mu_plus = TLorentzVector         (mu_plus_PX_DTF[0],            mu_plus_PY_DTF[0],             mu_plus_PZ_DTF[0],             mu_plus_PE_DTF[0]);
        if(Kst2Kspiplus){
            /// Ks, pi+
            LorVec_DTF_pi_plus = TLorentzVector         (pi_plus_PX_DTF[0],            pi_plus_PY_DTF[0],             pi_plus_PZ_DTF[0],             pi_plus_PE_DTF[0]);
            LorVec_DTF_KS_pi_plus = TLorentzVector      (Ks_pi_plus_PX_DTF[0],         Ks_pi_plus_PY_DTF[0],          Ks_pi_plus_PZ_DTF[0],          Ks_pi_plus_PE_DTF[0]);
            LorVec_DTF_KS_pi_minus = TLorentzVector     (Ks_pi_minus_PX_DTF[0],        Ks_pi_minus_PY_DTF[0],         Ks_pi_minus_PZ_DTF[0],         Ks_pi_minus_PE_DTF[0]);
            LorVec_DTF_pi_plus_MisIDed.SetXYZM    (pi_plus_PX_DTF[0],            pi_plus_PY_DTF[0],             pi_plus_PZ_DTF[0],             PDGMASS.MU);
            LorVec_DTF_mu_plus_MisIDed.SetXYZM    (mu_plus_PX_DTF[0],            mu_plus_PY_DTF[0],             mu_plus_PZ_DTF[0],             PDGMASS.PI_PLUS);
        }
        else{
            /// K+, pi0
            LorVec_DTF_K_plus = TLorentzVector    (K_plus_PX_DTF[0],             K_plus_PY_DTF[0],              K_plus_PZ_DTF[0],              K_plus_PE_DTF[0]);
            LorVec_DTF_gamma1 = TLorentzVector    (gamma1_PX_DTF[0],             gamma1_PY_DTF[0],              gamma1_PZ_DTF[0],              gamma1_PE_DTF[0]);
            LorVec_DTF_gamma2 = TLorentzVector    (gamma2_PX_DTF[0],             gamma2_PY_DTF[0],              gamma2_PZ_DTF[0],              gamma2_PE_DTF[0]);
        }
        //directly recombine the non-stable particles for later use (K*+ not really needed/used for Kst2Kspiplus)
        if(Kst2Kspiplus){
            LorVec_DTF_K_short = LorVec_DTF_KS_pi_plus + LorVec_DTF_KS_pi_minus;
            LorVec_DTF_KaonPion = LorVec_DTF_K_short + LorVec_DTF_pi_plus;
        }
        if(Kst2Kpluspi0Resolved){
            LorVec_DTF_DiGamma = LorVec_DTF_gamma1 + LorVec_DTF_gamma2;
            LorVec_DTF_pi_zero_resolved = LorVec_DTF_DiGamma;
            LorVec_DTF_KaonPion = LorVec_DTF_K_plus + LorVec_DTF_pi_zero_resolved;
        }
	LorVec_DTF_K_star_plus = LorVec_DTF_KaonPion;
        LorVec_DTF_DiMuon = LorVec_DTF_mu_minus + LorVec_DTF_mu_plus;
        LorVec_DTF_B_plus = LorVec_DTF_K_star_plus + LorVec_DTF_DiMuon;


        /////////////////////
        //    PDG values   //
        /////////////////////
        LorVec_PDG_B_plus.SetXYZM           (B_plus_PX,             B_plus_PY,              B_plus_PZ,              PDGMASS.B_PLUS);
        LorVec_PDG_K_star_plus.SetXYZM      (K_star_plus_PX,        K_star_plus_PY,         K_star_plus_PZ,         PDGMASS.K_STAR_PLUS);
        LorVec_PDG_mu_minus.SetXYZM         (mu_minus_PX,           mu_minus_PY,            mu_minus_PZ,            PDGMASS.MU);
        LorVec_PDG_mu_plus.SetXYZM          (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             PDGMASS.MU);
        if(Kst2Kspiplus){
            /// Ks, pi+
            LorVec_PDG_K_short.SetXYZM          (K_short_PX,            K_short_PY,             K_short_PZ,             PDGMASS.K_SHORT);
            LorVec_PDG_pi_plus.SetXYZM          (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             PDGMASS.PI_PLUS);
            LorVec_PDG_KS_pi_plus.SetXYZM       (Ks_pi_plus_PX,         Ks_pi_plus_PY,          Ks_pi_plus_PZ,          PDGMASS.PI_PLUS);
            LorVec_PDG_KS_pi_minus.SetXYZM      (Ks_pi_minus_PX,        Ks_pi_minus_PY,         Ks_pi_minus_PZ,         PDGMASS.PI_PLUS);

            LorVec_PDG_pi_plus_MisIDed.SetXYZM  (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             PDGMASS.MU);
            LorVec_PDG_mu_plus_MisIDed.SetXYZM  (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             PDGMASS.PI_PLUS);
        }
        else{
            /// K+, pi0
            LorVec_PDG_K_plus.SetXYZM           (K_plus_PX,             K_plus_PY,              K_plus_PZ,              PDGMASS.K_PLUS);
            LorVec_PDG_pi_zero_merged.SetXYZM   (pi_zero_merged_PX,     pi_zero_merged_PY,      pi_zero_merged_PZ,      PDGMASS.PI_ZERO);
            LorVec_PDG_pi_zero_resolved.SetXYZM (pi_zero_resolved_PX,   pi_zero_resolved_PY,    pi_zero_resolved_PZ,    PDGMASS.PI_ZERO);
            LorVec_PDG_gamma1.SetXYZM           (gamma1_PX,             gamma1_PY,              gamma1_PZ,              PDGMASS.GAMMA);
            LorVec_PDG_gamma2.SetXYZM           (gamma2_PX,             gamma2_PY,              gamma2_PZ,              PDGMASS.GAMMA);
        }
        //Recombination of 'Parent'-Particles
        LorVec_PDG_DiMuon = LorVec_PDG_mu_minus + LorVec_PDG_mu_plus;
        if(Kst2Kspiplus){
            LorVec_PDG_KaonPion = LorVec_PDG_K_short + LorVec_PDG_pi_plus;
        }
        if(Kst2Kpluspi0Merged){
            LorVec_PDG_KaonPion = LorVec_PDG_K_plus + LorVec_PDG_pi_zero_merged;
        }
        if(Kst2Kpluspi0Resolved){ //todo: save both dtf and fixed masses
            LorVec_PDG_DiGamma = LorVec_PDG_gamma1 + LorVec_PDG_gamma2;
            LorVec_PDG_KaonPion = LorVec_PDG_K_plus + LorVec_PDG_DiGamma;
        }

        //////////////////////
        //  detector values //
        //////////////////////
        LorVec_B_plus = TLorentzVector        (B_plus_PX,             B_plus_PY,              B_plus_PZ,              B_plus_PE);
        LorVec_K_star_plus = TLorentzVector   (K_star_plus_PX,        K_star_plus_PY,         K_star_plus_PZ,         K_star_plus_PE);
        LorVec_mu_minus = TLorentzVector      (mu_minus_PX,           mu_minus_PY,            mu_minus_PZ,            mu_minus_PE);
        LorVec_mu_plus = TLorentzVector       (mu_plus_PX,            mu_plus_PY,             mu_plus_PZ,             mu_plus_PE);
        if(Kst2Kspiplus){
            /// Ks, pi+
            LorVec_K_short = TLorentzVector       (K_short_PX,            K_short_PY,             K_short_PZ,             K_short_PE);
            LorVec_pi_plus = TLorentzVector       (pi_plus_PX,            pi_plus_PY,             pi_plus_PZ,             pi_plus_PE);
            LorVec_KS_pi_plus = TLorentzVector    (Ks_pi_plus_PX,         Ks_pi_plus_PY,          Ks_pi_plus_PZ,          Ks_pi_plus_PE);
            LorVec_KS_pi_minus = TLorentzVector   (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 = TLorentzVector            (K_plus_PX,             K_plus_PY,              K_plus_PZ,              K_plus_PE);
            LorVec_pi_zero_merged = TLorentzVector	(pi_zero_merged_PX,     pi_zero_merged_PY,      pi_zero_merged_PZ,      pi_zero_merged_PE);
            LorVec_pi_zero_resolved = TLorentzVector  (pi_zero_resolved_PX,   pi_zero_resolved_PY,    pi_zero_resolved_PZ,    pi_zero_resolved_PE);
            LorVec_gamma1 = TLorentzVector            (gamma1_PX,             gamma1_PY,              gamma1_PZ,              gamma1_PE);
            LorVec_gamma2 = TLorentzVector            (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_DTF_PDG_KshortMuMu   = LorVec_DTF_PDG_K_short+ LorVec_DTF_PDG_DiMuon;
            LorVec_DTF_KshortMuMu       = LorVec_DTF_K_short    + LorVec_DTF_DiMuon;
            LorVec_PDG_KshortMuMu       = LorVec_PDG_K_short    + LorVec_PDG_DiMuon;
            LorVec_KshortMuMu           = LorVec_K_short        + LorVec_DiMuon;
        }
        if(Kst2Kpluspi0Merged){
            LorVec_DTF_PDG_KplusMuMu    = LorVec_DTF_PDG_K_plus + LorVec_DTF_PDG_DiMuon;
            LorVec_DTF_KplusMuMu        = LorVec_DTF_K_plus     + LorVec_DTF_DiMuon;
            LorVec_PDG_KplusMuMu        = LorVec_PDG_K_plus     + LorVec_PDG_DiMuon;
            LorVec_KplusMuMu            = LorVec_K_plus         + LorVec_DiMuon;
        }
        if(Kst2Kpluspi0Resolved){
            LorVec_DTF_PDG_KplusMuMu    = LorVec_DTF_PDG_K_plus + LorVec_DTF_PDG_DiMuon;
            LorVec_DTF_KplusMuMu        = LorVec_DTF_K_plus     + LorVec_DTF_DiMuon;
            LorVec_PDG_KplusMuMu        = LorVec_PDG_K_plus     + LorVec_PDG_DiMuon;
            LorVec_KplusMuMu            = LorVec_K_plus         + LorVec_DiMuon;
            LorVec_DTF_PDG_ReconBPiZero = LorVec_DTF_PDG_DiMuon + LorVec_DTF_PDG_K_plus + LorVec_DTF_PDG_pi_zero_resolved;
            LorVec_DTF_ReconBPiZero     = LorVec_DTF_DiMuon     + LorVec_DTF_K_plus     + LorVec_DTF_pi_zero_resolved;
            LorVec_PDG_ReconBPiZero     = LorVec_PDG_DiMuon     + LorVec_PDG_K_plus     + LorVec_PDG_pi_zero_resolved;
            LorVec_ReconBPiZero         = LorVec_DiMuon         + LorVec_K_plus         + LorVec_pi_zero_resolved;
            LorVec_DTF_PDG_pi_zero_resolved  = LorVec_DTF_PDG_DiGamma;
            LorVec_DTF_pi_zero_resolved      = LorVec_DTF_DiGamma;
            LorVec_PDG_pi_zero_resolved      = LorVec_PDG_DiGamma;
            LorVec_pi_zero_resolved          = LorVec_DiGamma;
        }

        LorVec_DTF_PDG_K_star_plus  = LorVec_DTF_PDG_KaonPion;
        LorVec_DTF_K_star_plus      = LorVec_DTF_KaonPion;
        LorVec_PDG_K_star_plus      = LorVec_PDG_KaonPion;
        LorVec_K_star_plus          = LorVec_KaonPion;
        LorVec_DTF_PDG_ReconB   = LorVec_DTF_PDG_DiMuon + LorVec_DTF_PDG_K_star_plus;
        LorVec_DTF_ReconB       = LorVec_DTF_DiMuon     + LorVec_DTF_K_star_plus;
        LorVec_PDG_ReconB       = LorVec_PDG_DiMuon     + LorVec_PDG_K_star_plus;
        LorVec_ReconB           = LorVec_DiMuon         + LorVec_K_star_plus;


        //fill MASS HISTOGRAMS
        h_B_M_DTF_PDG   ->Fill(LorVec_DTF_PDG_ReconB.M());
        h_B_M_DTF       ->Fill(LorVec_DTF_ReconB.M());
        h_B_M_PDG       ->Fill(LorVec_PDG_ReconB.M());
        h_B_M           ->Fill(LorVec_ReconB.M());
        h_Kstar_M_DTF_PDG   ->Fill(LorVec_DTF_PDG_K_star_plus.M());
        h_Kstar_M_DTF       ->Fill(LorVec_DTF_K_star_plus.M());
        h_Kstar_M_PDG       ->Fill(LorVec_PDG_K_star_plus.M());
        h_Kstar_M           ->Fill(LorVec_K_star_plus.M());
        if(Kst2Kspiplus){
            h_Kshort_M_DTF_PDG   ->Fill(LorVec_DTF_PDG_K_short.M());
            h_Kshort_M_DTF       ->Fill(LorVec_DTF_K_short.M());
            h_Kshort_M_PDG       ->Fill(LorVec_PDG_K_short.M());
            h_Kshort_M           ->Fill(LorVec_K_short.M());
        }
        if (Kst2Kpluspi0Resolved){
            h_pi0_M_DTF_PDG ->Fill(LorVec_DTF_PDG_pi_zero_resolved.M());
            h_pi0_M_DTF     ->Fill(LorVec_DTF_pi_zero_resolved.M());
            h_pi0_M_PDG     ->Fill(LorVec_PDG_pi_zero_resolved.M());
            h_pi0_M         ->Fill(LorVec_pi_zero_resolved.M());
        }

        //fill PT HISTOGRAMS
        h_B_PT_DTF_PDG   ->Fill(LorVec_DTF_PDG_ReconB.Pt());
        h_B_PT_DTF       ->Fill(LorVec_DTF_ReconB.Pt());
        h_B_PT_PDG       ->Fill(LorVec_PDG_ReconB.Pt());
        h_B_PT           ->Fill(LorVec_ReconB.Pt());
        h_Kstar_PT_DTF_PDG  ->Fill(LorVec_DTF_PDG_K_star_plus.Pt());
        h_Kstar_PT_DTF      ->Fill(LorVec_DTF_K_star_plus.Pt());
        h_Kstar_PT_PDG      ->Fill(LorVec_PDG_K_star_plus.Pt());
        h_Kstar_PT          ->Fill(LorVec_K_star_plus.Pt());
        if(Kst2Kspiplus){
            h_Kshort_PT_DTF_PDG  ->Fill(LorVec_DTF_PDG_K_short.Pt());
            h_Kshort_PT_DTF      ->Fill(LorVec_DTF_K_short.Pt());
            h_Kshort_PT_PDG      ->Fill(LorVec_PDG_K_short.Pt());
            h_Kshort_PT          ->Fill(LorVec_K_short.Pt());
        }
        if (Kst2Kpluspi0Resolved){
            h_pi0_PT_DTF_PDG ->Fill(LorVec_DTF_PDG_pi_zero_resolved.Pt());
            h_pi0_PT_DTF     ->Fill(LorVec_DTF_pi_zero_resolved.Pt());
            h_pi0_PT_PDG     ->Fill(LorVec_PDG_pi_zero_resolved.Pt());
            h_pi0_PT         ->Fill(LorVec_pi_zero_resolved.Pt());
        }
        //fill ETA HISTOGRAMS
        h_B_ETA_DTF_PDG   ->Fill(LorVec_DTF_PDG_ReconB.Eta() );
        h_B_ETA_DTF       ->Fill(LorVec_DTF_ReconB.Eta());
        h_B_ETA_PDG       ->Fill(LorVec_PDG_ReconB.Eta());
        h_B_ETA           ->Fill(LorVec_ReconB.Eta());
        h_Kstar_ETA_DTF_PDG   ->Fill(LorVec_DTF_PDG_K_star_plus.Eta());
        h_Kstar_ETA_DTF       ->Fill(LorVec_DTF_K_star_plus.Eta());
        h_Kstar_ETA_PDG       ->Fill(LorVec_PDG_K_star_plus.Eta());
        h_Kstar_ETA           ->Fill(LorVec_K_star_plus.Eta());
        if(Kst2Kspiplus){
            h_Kshort_ETA_DTF_PDG   ->Fill(LorVec_DTF_PDG_K_short.Eta());
            h_Kshort_ETA_DTF       ->Fill(LorVec_DTF_K_short.Eta());
            h_Kshort_ETA_PDG       ->Fill(LorVec_PDG_K_short.Eta());
            h_Kshort_ETA           ->Fill(LorVec_K_short.Eta());
        }
        if (Kst2Kpluspi0Resolved){
            h_pi0_ETA_DTF_PDG ->Fill(LorVec_DTF_PDG_pi_zero_resolved.Eta());
            h_pi0_ETA_DTF     ->Fill(LorVec_DTF_pi_zero_resolved.Eta());
            h_pi0_ETA_PDG     ->Fill(LorVec_PDG_pi_zero_resolved.Eta());
            h_pi0_ETA         ->Fill(LorVec_pi_zero_resolved.Eta());
        }
        //fill angles
        h_ctk->Fill(ThetaK);
        h_ctl->Fill(ThetaL);
        h_sp->Fill(Phi);


        //fill other variables
        if (Kst2Kpluspi0Resolved){
            OpenAngle_DTF_PDG  = LorVec_DTF_PDG_gamma1.Angle(LorVec_DTF_PDG_gamma2.Vect());
            OpenAngle_DTF      = LorVec_DTF_gamma1.Angle(LorVec_DTF_gamma2.Vect());
            OpenAngle_PDG      = LorVec_PDG_gamma1.Angle(LorVec_PDG_gamma2.Vect());
            OpenAngle          = LorVec_gamma1.Angle(LorVec_gamma2.Vect());
            h_OpenAngle_DTF_PDG ->Fill(OpenAngle_DTF_PDG);
            h_OpenAngle_DTF     ->Fill(OpenAngle_DTF);
            h_OpenAngle_PDG     ->Fill(OpenAngle_PDG);
            h_OpenAngle         ->Fill(OpenAngle);

        }


    } //end of event loop

    // Kst2Kpluspi0Resolved case
    string newName = "";
    preselected ?  newName = "_DTFcheck_preselected" : newName ="_DTFcheck";
    if (truthMatched) newName = newName + "_TM";

    if(Kst2Kpluspi0Resolved){
        if(!MC){
            output = new TFile(Form("%s/data/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
        }
        else{
            if(ReferenceChannel){
                output = new TFile(Form("%s/data/MC/RefKplusPi0/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
            else if(PHSP){
                output = new TFile(Form("%s/data/PHSP/KplusPi0/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
            else{
                output = new TFile(Form("%s/data/MC/KplusPi0/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
        }
    }
    // Kst2Kpluspi0Merged case
    if(Kst2Kpluspi0Merged){
        if(!MC){
            output = new TFile(Form("%s/data/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
        }
        else{
            if(ReferenceChannel){
                output = new TFile(Form("%s/data/MC/RefKplusPi0/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
            else if(PHSP){
                output = new TFile(Form("%s/data/PHSP/KplusPi0/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
            else{
                output = new TFile(Form("%s/data/MC/KplusPi0/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
        }
    }

    // Kst2Kspiplus case
    if(Kst2Kspiplus){
        if(!MC){
            output = new TFile(Form("%s/data/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
        }
        else{
            if(ReferenceChannel){
                output = new TFile(Form("%s/data/MC/RefKshortPiplus/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
            else if(PHSP){
                output = new TFile(Form("%s/data/PHSP/KshortPiplus/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
            else{
                output = new TFile(Form("%s/data/MC/KshortPiplus/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()),"RECREATE");
            }
        }
    }

    output->cd();

    //save MASS HISTOGRAMS
    h_B_M_DTF_PDG   ->Write();
    h_B_M_DTF       ->Write();
    h_B_M_PDG       ->Write();
    h_B_M           ->Write();
    h_Kstar_M_DTF_PDG   ->Write();
    h_Kstar_M_DTF       ->Write();
    h_Kstar_M_PDG       ->Write();
    h_Kstar_M           ->Write();
    if(Kst2Kspiplus){
        h_Kshort_M_DTF_PDG   ->Write();
        h_Kshort_M_DTF       ->Write();
        h_Kshort_M_PDG       ->Write();
        h_Kshort_M           ->Write();
    }
    if (Kst2Kpluspi0Resolved){
        h_pi0_M_DTF_PDG     ->Write();
        h_pi0_M_DTF         ->Write();
        h_pi0_M_PDG         ->Write();
        h_pi0_M             ->Write();
    }

    //save PT HISTOGRAMS
    h_B_PT_DTF_PDG  ->Write();
    h_B_PT_DTF      ->Write();
    h_B_PT_PDG      ->Write();
    h_B_PT          ->Write();
    h_Kstar_PT_DTF_PDG  ->Write();
    h_Kstar_PT_DTF      ->Write();
    h_Kstar_PT_PDG      ->Write();
    h_Kstar_PT          ->Write();
    if(Kst2Kspiplus){
        h_Kshort_PT_DTF_PDG  ->Write();
        h_Kshort_PT_DTF      ->Write();
        h_Kshort_PT_PDG      ->Write();
        h_Kshort_PT          ->Write();
    }
    if (Kst2Kpluspi0Resolved){
        h_pi0_PT_DTF_PDG    ->Write();
        h_pi0_PT_DTF        ->Write();
        h_pi0_PT_PDG        ->Write();
        h_pi0_PT            ->Write();
    }

    //save ETA HISTOGRAMS
    h_B_ETA_DTF_PDG ->Write();
    h_B_ETA_DTF     ->Write();
    h_B_ETA_PDG     ->Write();
    h_B_ETA         ->Write();
    h_Kstar_ETA_DTF_PDG ->Write();
    h_Kstar_ETA_DTF     ->Write();
    h_Kstar_ETA_PDG     ->Write();
    h_Kstar_ETA         ->Write();
    if(Kst2Kspiplus){
        h_Kshort_ETA_DTF_PDG ->Write();
        h_Kshort_ETA_DTF     ->Write();
        h_Kshort_ETA_PDG     ->Write();
        h_Kshort_ETA         ->Write();
    }
    if (Kst2Kpluspi0Resolved){
        h_pi0_ETA_DTF_PDG   ->Write();
        h_pi0_ETA_DTF       ->Write();
        h_pi0_ETA_PDG       ->Write();
        h_pi0_ETA           ->Write();

        //save other histograms
        h_OpenAngle_DTF_PDG ->Write();
        h_OpenAngle_DTF     ->Write();
        h_OpenAngle_PDG     ->Write();
        h_OpenAngle         ->Write();
    }

    //save Angle histograms
//    h_ctk   ->Write();
//    h_ctl   ->Write();
//    h_sp    ->Write();

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

int DTFcomparisonAllData(Int_t Run = 1, bool preselected = true) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }

    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"down", false, false, false, preselected, 0) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , false, false, false, preselected, 0) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"down", false, false, false, preselected, 0) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , false, false, false, preselected, 0) == 0) return 0;
    return 1;

}

int DTFcomparisonAllMC(Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }

    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"down", true, false, false, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , true, false, false, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"down", true, false, false, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , true, false, false, preselected, truthMatched) == 0) return 0;
    return 1;

}

int DTFcomparisonAllRefChannel(Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }
    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"down", true, true, false, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , true, true, false, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"down", true, true, false, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , true, true, false, preselected, truthMatched) == 0) return 0;
    return 1;

}

int DTFcomparisonAllPHSP(Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }
    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"down", true, false, true, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , true, false, true, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"down", true, false, true, preselected, truthMatched) == 0) return 0;
    if ( DTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , true, false, true, preselected, truthMatched) == 0) return 0;
    return 1;

}

int DTFcomparisonAll( Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if (DTFcomparisonAllData    (Run, preselected)==0) return 0;
    if (DTFcomparisonAllMC      (Run, preselected, truthMatched)==0) return 0;
    if (Run == 2) cout << "[INFO]\tNo reference channel for run 2, Skipping.";
    else{
         if (DTFcomparisonAllRefChannel  (Run, preselected, truthMatched)==0) return 0;
    }
    if (DTFcomparisonAllPHSP    (Run, preselected, truthMatched)==0) return 0;
    return 1;
}

int DTFcomparisonComplete( Int_t Run = 1){
    if (DTFcomparisonAllData    (Run, 1)==0) return 0;
    if (DTFcomparisonAllData    (Run, 0)==0) return 0;
    if (DTFcomparisonAllMC          (Run, 0, 0)==0) return 0;
    if (DTFcomparisonAllMC          (Run, 1, 0)==0) return 0;
    if (DTFcomparisonAllMC          (Run, 1, 1)==0) return 0;
    if (DTFcomparisonAllRefChannel  (Run, 0, 0)==0) return 0;
    if (DTFcomparisonAllRefChannel  (Run, 1, 0)==0) return 0;
    if (DTFcomparisonAllRefChannel  (Run, 1, 1)==0) return 0;
    if (DTFcomparisonAllPHSP        (Run, 0, 0)==0) return 0;
    if (DTFcomparisonAllPHSP        (Run, 1, 0)==0) return 0;
    if (DTFcomparisonAllPHSP        (Run, 1, 1)==0) return 0;
    return 1;
}


int plotDTFcomparison(std::string year = "2011", std::string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false, bool preselected = true, bool truthMatched = false){

    // Kst2Kpluspi0Resolved case
    TFile *input = 0;
    string newName = "";
    string legTitle = "";
    string TM = "";
    preselected ?  newName = "_DTFcheck_preselected" : newName ="_DTFcheck";
    truthMatched ?  TM = "_TM" : "";
    newName = newName + TM;
    Kst2Kspiplus?  legTitle= particleLatex("Kstar") +" #rightarrow " +particleLatex("Kshort") + particleLatex("piplus") :legTitle= particleLatex("Kstar") +" #rightarrow " +particleLatex("Kplus") + particleLatex("pizero");
    preselected ?  legTitle= legTitle + ": preselected" : legTitle = legTitle+": stripped";
    truthMatched ?  legTitle= "#splitline{" + legTitle+":}{TruthMatched}" : legTitle = legTitle;
    if (MC){
        if(ReferenceChannel) legTitle = legTitle +  " {MCRef}";
        else if (PHSP) legTitle = legTitle +  " {PHSP MC}";
        else legTitle = legTitle +  " {MC}";
    }


    if(Kst2Kpluspi0Resolved){
        if(!MC){
            input = TFile::Open(Form("%s/data/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            cout << Form("%s/data/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()) << endl;
        }
        else{
            if(ReferenceChannel){
                input = TFile::Open(Form("%s/data/MC/RefKplusPi0/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
            else if(PHSP){
                input = TFile::Open(Form("%s/data/PHSP/KplusPi0/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
            else{
                input = TFile::Open(Form("%s/data/MC/KplusPi0/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
        }
    }
    // Kst2Kpluspi0Merged case
    if(Kst2Kpluspi0Merged){
        if(!MC){
            input = TFile::Open(Form("%s/data/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
        }
        else{
            if(ReferenceChannel){
                input = TFile::Open(Form("%s/data/MC/RefKplusPi0/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
            else if(PHSP){
                input = TFile::Open(Form("%s/data/PHSP/KplusPi0/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
            else{
                input = TFile::Open(Form("%s/data/MC/KplusPi0/%s%s/%s%s_pi0Merged%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
        }
    }
    // Kst2Kspiplus case
    if(Kst2Kspiplus){
        if(!MC){
            input = TFile::Open(Form("%s/data/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
        }
        else{
            if(ReferenceChannel){
                input = TFile::Open(Form("%s/data/MC/RefKshortPiplus/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
            else if(PHSP){
                input = TFile::Open(Form("%s/data/PHSP/KshortPiplus/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
            else{
                input = TFile::Open(Form("%s/data/MC/KshortPiplus/%s%s/%s%s_piplus%s.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),newName.c_str()));
            }
        }
    }


    string tmpName = "";
    string tmpNameCanvas = "";
    int col = 1;
    if(Kst2Kpluspi0Resolved){
        for (auto& var: Variables){
            for (auto& particle: ParticlesPiZero){
                tmpNameCanvas = "c_" + particle + "_" + var;
                TCanvas *canvas = new TCanvas(tmpNameCanvas.c_str(), "", 10,10,600,600);
                designCanvas(canvas);
                canvas->cd();
                gStyle->SetOptStat(0);
                TLegend *leg = new TLegend(0.57,0.78,0.96,0.96);
                col = 1;
                for (auto& opt: DTFoptions){
                    tmpName = "h_" + particle + "_" + var + "_" + opt;
                    if (strcmp(&tmpName[tmpName.length()-1],"_")==0) tmpName.pop_back(); //if ending with _, remove it
                    TH1D* tmp =  (TH1D*)input->Get(tmpName.c_str());  //check memory leaks
                    designPlots(tmp,col);
                    designLegend(leg, legTitle, tmp, opt);
                    //tmp->SetTitle((variableLatex(var)+"(" + particleLatex(particle)+")").c_str());
                    tmp->Draw("SAME");
                    col++;
                }
                leg->Draw();
                tmpName = path_to_output_KplusPizero;
                if (MC){
                    if(ReferenceChannel) tmpName = tmpName + "/data/MC/RefKplusPi0/" + year + magnet;
                    else if (PHSP) tmpName  = tmpName + "/data/PHSP/KplusPi0/" + year + magnet;
                    else tmpName  = tmpName + "/data/MC/KplusPi0/" + year + magnet;
                }
                else tmpName  = tmpName + "/data/" + year + magnet;
                tmpNameCanvas = tmpName + "/" + tmpNameCanvas + TM +  (preselected ? "_preselected" :"") + "."  +figureFormat;
                canvas->Print(tmpNameCanvas.c_str(),figureFormat.c_str());             //#TODO: cancel piplus, replace by path_to accordingly
                leg->Clear();
                canvas->Clear();
                delete(leg);
                delete(canvas);
            }
        }
        for (auto& evar: ExtraVariables){
            tmpNameCanvas = "c_" + evar;
            TCanvas *canvas = new TCanvas(tmpNameCanvas.c_str(), "", 10,10,600,600);
            designCanvas(canvas);
            canvas->cd();
            gStyle->SetOptStat(0);
            TLegend *leg = new TLegend(0.62,0.8,0.96,0.96);
            col = 1;
            for (auto& opt: DTFoptions){
                tmpName = "h_" + evar + "_" + opt;
                if (strcmp(&tmpName[tmpName.length()-1],"_")==0) tmpName.pop_back(); //if ending with _, remove it
                TH1D* tmp =  (TH1D*)input->Get(tmpName.c_str());  //check memory leaks
                designPlots(tmp,col);
                designLegend(leg,  legTitle, tmp, opt);
               // tmp->SetTitle((variableLatex(evar)).c_str());
                tmp->Draw("SAME");
                col++;
            }
            leg->Draw();
            tmpName = path_to_output_KplusPizero;
            if (MC){
                if(ReferenceChannel) tmpName = tmpName + "/data/MC/RefKplusPi0/" + year + magnet;
                else if (PHSP) tmpName  = tmpName + "/data/PHSP/KplusPi0/" + year + magnet;
                else tmpName  = tmpName + "/data/MC/KplusPi0/" + year + magnet;
            }
            else tmpName  = tmpName + "/data/" + year + magnet;
            tmpNameCanvas = tmpName + "/" + tmpNameCanvas + TM +  (preselected ? "_preselected" :"") + "."  +figureFormat;
            canvas->Print(tmpNameCanvas.c_str(),figureFormat.c_str());             //#TODO: cancel piplus, replace by path_to accordingly
            leg->Clear();
            canvas->Clear();
            delete(leg);
            delete(canvas);

        }
    }
    if(Kst2Kspiplus){
        for (auto& var: Variables){
            for (auto& particle: ParticlesKShort){
                tmpNameCanvas = "c_" + particle + "_" + var;
                TCanvas *canvas = new TCanvas(tmpNameCanvas.c_str(), "", 10,10,600,600);
                designCanvas(canvas);
                canvas->cd();
                gStyle->SetOptStat(0);
                TLegend *leg = new TLegend(0.62,0.8,0.96,0.96);
                col = 1;
                for (auto& opt: DTFoptions){
                    tmpName = "h_" + particle + "_" + var + "_" + opt;
                    if (strcmp(&tmpName[tmpName.length()-1],"_")==0) tmpName.pop_back(); //if ending with _, remove it
                    TH1D* tmp =  (TH1D*)input->Get(tmpName.c_str());  //check memory leaks
                    designPlots(tmp,col);
                    designLegend(leg,  legTitle, tmp, opt);
                  //  tmp->SetTitle((variableLatex(var)+"(" + particleLatex(particle)+")").c_str());
                    tmp->Draw("SAME");
                    col++;
                }
                leg->Draw();
                tmpName = path_to_output_KshortPiplus;
                if (MC){
                    if(ReferenceChannel) tmpName = tmpName + "/data/MC/RefKplusPi0/" + year + magnet;
                    else if (PHSP) tmpName  = tmpName + "/data/PHSP/KplusPi0/" + year + magnet;
                    else tmpName  = tmpName + "/data/MC/KplusPi0/" + year + magnet;
                }
                else tmpName  = tmpName + "/data/" + year + magnet;
                tmpNameCanvas = tmpName + "/piplus_" + tmpNameCanvas +TM +  (preselected ? "_preselected" :"") +  "." +figureFormat;
                canvas->Print(tmpNameCanvas.c_str(),figureFormat.c_str());             //#TODO: cancel piplus, replace by path_to accordingly
                leg->Clear();
                canvas->Clear();
                delete(leg);
                delete(canvas);
            }
        }
    }
    return 1;
}

int plotDTFcomparisonAllData(Int_t Run = 1, bool preselected = true) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }

    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"down", false, false, false, preselected, 0) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , false, false, false, preselected, 0) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"down", false, false, false, preselected, 0) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , false, false, false, preselected, 0) == 0) return 0;
    return 1;

}

int plotDTFcomparisonAllMC(Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }

    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"down", true, false, false, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , true, false, false, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"down", true, false, false, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , true, false, false, preselected, truthMatched) == 0) return 0;
    return 1;

}

int plotDTFcomparisonAllRefChannel(Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }
    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"down", true, true, false, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , true, true, false, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"down", true, true, false, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , true, true, false, preselected, truthMatched) == 0) return 0;
    return 1;

}

int plotDTFcomparisonAllPHSP(Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if(Run != 1 && Run != 2){
        std::cout << "[ERROR]\t\tInvalid Run number given: " << Run  << ". Exit program!" << std::endl;
        return 0;
    }
    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"down", true, false, true, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2011" : "2015"),"up"  , true, false, true, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"down", true, false, true, preselected, truthMatched) == 0) return 0;
    if ( plotDTFcomparison((Run == 1 ? "2012" : "2016"),"up"  , true, false, true, preselected, truthMatched) == 0) return 0;
    return 1;

}

int plotDTFcomparisonAll( Int_t Run = 1, bool preselected = true, bool truthMatched = false) {

    if (plotDTFcomparisonAllData        (Run, preselected)==0) return 0;
    if (plotDTFcomparisonAllMC          (Run, preselected, truthMatched)==0) return 0;
    if (plotDTFcomparisonAllRefChannel  (Run, preselected, truthMatched)==0) return 0;
    if (plotDTFcomparisonAllPHSP        (Run, preselected, truthMatched)==0) return 0;
    return 1;
}

int plotDTFcomparisonComplete( Int_t Run = 1){
    if (plotDTFcomparisonAllData    (Run, 1)==0) return 0;
    if (plotDTFcomparisonAllData    (Run, 0)==0) return 0;
    if (plotDTFcomparisonAllMC          (Run, 0, 0)==0) return 0;
    if (plotDTFcomparisonAllMC          (Run, 1, 0)==0) return 0;
    if (plotDTFcomparisonAllMC          (Run, 1, 1)==0) return 0;
    if (plotDTFcomparisonAllRefChannel  (Run, 0, 0)==0) return 0;
    if (plotDTFcomparisonAllRefChannel  (Run, 1, 0)==0) return 0;
    if (plotDTFcomparisonAllRefChannel  (Run, 1, 1)==0) return 0;
    if (plotDTFcomparisonAllPHSP        (Run, 0, 0)==0) return 0;
    if (plotDTFcomparisonAllPHSP        (Run, 1, 0)==0) return 0;
    if (plotDTFcomparisonAllPHSP        (Run, 1, 1)==0) return 0;
    return 1;
}
