EWP-BplusToKstMuMu-AngAna/Code/Selection/BDTcutScanner.cpp

//Functions that calculate/fit the yield for given MVA response cut
//and make pretty plots from it
//Renata Kopecna

#include "GlobalFunctions.hh"
#include "Paths.hpp"
#include "MassFit.hpp"
#include "Utils.hpp"
#include "./BmassShape/SignalType.hpp"
#include "EfficiencyClass.cpp"


using namespace std;
using namespace RooFit;
using namespace RooStats;

class YieldInfo{
 public:
    string year;
    double sigYield;
    double sigYieldErr;
    double refYield;
    double refYieldErr;
    double bkgYield;
    double bkgYieldErr;
    double sigEff;
    double refEff;
    double allSigEvts;
    YieldInfo(){ //default constructor
        year        = "";
        sigYield    = 0;
        sigYieldErr = 0;
        refYield    = 0;
        refYieldErr = 0;
        bkgYield    = 0;
        bkgYieldErr = 0;
        sigEff      = 0;
        refEff      = 0;
        allSigEvts  = 0;
    }
    void addYield(YieldInfo addInfo){
        sigYield    += addInfo.sigYield;
        sigYieldErr += addInfo.sigYieldErr;
        refYield    += addInfo.refYield;
        refYieldErr += addInfo.refYieldErr;
        bkgYield    += addInfo.bkgYield;
        bkgYieldErr += addInfo.bkgYieldErr;
        sigEff      += addInfo.sigEff; //this doesn't make much sense
        refEff      += addInfo.refEff; //this doesn't make much sense
        allSigEvts  += addInfo.allSigEvts;
    }

    ~YieldInfo();   //destuctor
};

YieldInfo::~YieldInfo(){//destuctor

}

bool FixShape = true;

bool RemoveMultiple = true; //For now use as a global boolean

//Function that takes upper mass sideband histogram and fits it with an exponential
//Using a Tree, starting the upper mass sideband from range+PDGMASS.Bplus and possibly cutting away stuff defined by cut
double GetBackgroundFromSidebandFit(TChain *Tree, double range, string cut, string path){

    TString massBranch = UseDTF ? "B_plus_M" : "B_plus_DTF_M";
    bool savePlots = true;

    //Fill histogram from Tree
    int nBins = 50;
    TH1D *h_B_mass = new TH1D("h_B_mass","h_B_mass",nBins,PDGMASS.B_PLUS + range, cut_B_plus_M_high);
    double binWidth = (cut_B_plus_M_high - PDGMASS.B_PLUS - range)/nBins;
    Tree->Draw(massBranch+" >> h_B_mass",cut.c_str());

    //Define fit function
    TF1 *f_bkg = new TF1("f_bkg","expo",cut_B_plus_M_low,cut_B_plus_M_high);
    design_function(f_bkg,kRed,1);

    //Fit
    h_B_mass->Fit(f_bkg,"Q","", PDGMASS.B_PLUS + range,cut_B_plus_M_high); //exponential is f(x) = exp(p0+p1*x)
    double integral = f_bkg->Integral(PDGMASS.B_PLUS-range,PDGMASS.B_PLUS+range)/binWidth;

    //Save the plots
    if (savePlots){
        string name = path.substr(path.find("/Kplus")+6,path.length());
        //Create canvas
        TCanvas  *c = c_canvas(name);
        c->cd();
        design_TH1D(h_B_mass,"B^{+} mass [MeV]","Counts a.u.", kBlack);
        h_B_mass->GetYaxis()->SetTitleOffset(1.05);
        h_B_mass->Draw();
        //Create TPaveText
        TPaveText *fitParams= new TPaveText(0.6,0.7,0.87,0.9,"NDC");
        fitParams->SetFillColor(kWhite);
        fitParams->AddText(Form("Offset : %.2f ",  f_bkg->GetParameter(0)));
        fitParams->AddText(Form("Exp: %.2f%%", f_bkg->GetParameter(1)*100));
        fitParams->AddText(Form("Integral: %.2f ", integral));
        fitParams->Draw("SAME");

        //Create file
        TFile *fitFile  = new TFile(path.c_str(),"RECREATE");
        coutDebug("Opening " + string(fitFile->GetPath()));
        fitFile->cd();
        h_B_mass->Write();
        f_bkg->Write();
        c->Write();
        replace(path,".root",".eps");
        c->SaveAs(path.c_str(), ".eps");
        fitFile->Close();
    }
    //return the integral
    h_B_mass->Clear();
    return integral;

}

double GetMVAefficiency(string year, int Run, bool KshortDecaysInVelo,
                        bool UseLowQ2Range, bool IncludeMultipleEff, bool sig, double TMVAcut){
    //Load the efficiencies
    coutDebug("Openning MVA efficiency files.");
    string path = GetEfficiencyFile("BDT", year, "both", Run, !sig, sig,  false, KshortDecaysInVelo, IncludeMultipleEff, true, UseLowQ2Range, false,"");

    if (Run == 0 && std::stoi(year) > 2016 && !sig) replace(path,year,"2016"); //RefChann available for 2015-2016
    if (Run == 0 && year=="2015" && sig) replace(path,year,"2016");            //No reasonable MC for 2015

    coutDebug("Openning file: " + path);
    TFile * effFile = new TFile(path.c_str(), "OPEN");

    //Get the TGraphs
    coutDebug("File opened, getting the graphs now.");
    string graphName = "effScan_" + (Run ==0 ? year : "Run" + to_string(Run));
    if (Run == 0 && std::stoi(year) > 2016) replace(graphName,year,"2016");
    if (Run == 0 && year=="2015" && sig)    replace(graphName,year,"2016");
    coutDebug("\t" + graphName);

    TGraphErrors *effGraph= (TGraphErrors*)effFile->Get(TMVAmethod+graphName.c_str());
    double efficiency = effGraph->Eval(TMVAcut);

    coutDebug(string(sig ? "Signal   " : "Reference") +" MVA efficiency: " + to_string(efficiency));

    //Close the efficiency files
    effFile->Close();
    return efficiency;

}

double GetSelectionEfficiency(string year, int Run, bool KshortDecaysInVelo,
                              bool UseLowQ2Range, bool RemoveMultiple, bool sig){
    //Load the efficiencies
    coutDebug("Openning efficiency files.");
    EffAndError eff = getSelectionEfficiencySimple(true,year,Run,sig,!sig, false, KshortDecaysInVelo, RemoveMultiple, true, "", -1, "", gammaTMdefault);

    coutDebug(string(sig ? "Signal   " : "Reference") +" selection efficiency: " + to_string(eff.value));

    //Close the efficiency files
    return eff.value;
}


string GetBackgroundFunction(bool KshortDecaysInVelo){    //Set signal and background function separately for each decay channel
    string BackGroundFunction = "SingleExponential";
    if(Kst2Kspiplus){
        if(KshortDecaysInVelo) BackGroundFunction = "SingleExponential";
        else                   BackGroundFunction = "DoubleExponential";
    }
    return BackGroundFunction;
}

//Takes yield from the reference channel,
YieldInfo GetSigAndBkgEstimation(string year, int Run, bool KshortDecaysInVelo,
                                 bool UseLowQ2Range, Double_t TMVAcut, bool scan){

    bool fixedMassRegion = !Kst2Kspiplus; //Calculates the signal from +- 100 MeV or +- SignalRegionNsigma*effective sigma

    //initialize the struct
    YieldInfo yield = YieldInfo();
    yield.year = year;

    bool RemoveMultiple = true;

    //Set signal and background function separately for each decay channel
    string SignalFunction = "OneCB";
    string BackGroundFunction = GetBackgroundFunction(KshortDecaysInVelo);

    //Get the reference signal yield
    if (!scan){
        coutInfo("Running the fit!");
        yield.refYield = basicYieldFit(year,Run,
                                       false, false, true, false,
                                       FixShape, SignalFunction, BackGroundFunction, FixShape,
                                       KshortDecaysInVelo, UseLowQ2Range,
                                       TMVAcut, fixedMassRegion, true, RemoveMultiple);
    }

    //Load the file from fitting
    TFile *fitFile  = new TFile(GetMassFitFile(year, Run,
                                               false, false, true, false,
                                               FixShape, SignalFunction, BackGroundFunction, FixShape,
                                               KshortDecaysInVelo, UseLowQ2Range,
                                               TMVAcut, fixedMassRegion, RemoveMultiple).c_str(),"OPEN");

    coutDebug("Opening " + string(fitFile->GetPath()));

    if (scan) yield.refYield = getSigYield(fitFile);
    double effSigma = fixedMassRegion ? B_plus_M_signal_window : getEffSigma(fitFile);

    if (yield.refYield == 0) coutWarning("Reference channel yield is zero!");

    yield.sigEff = GetMVAefficiency(year,Run,KshortDecaysInVelo,UseLowQ2Range,RemoveMultiple,true,TMVAcut) * GetSelectionEfficiency(year, Run, KshortDecaysInVelo, UseLowQ2Range, !RemoveMultiple, true);
    yield.refEff = GetMVAefficiency(year,Run,KshortDecaysInVelo,UseLowQ2Range,RemoveMultiple,false,TMVAcut) * GetSelectionEfficiency(year, Run, KshortDecaysInVelo, UseLowQ2Range, !RemoveMultiple, false);

    coutDebug("Sig Yield efficiency = " + to_string(yield.sigEff));
    coutDebug("Ref Yield efficiency = " + to_string(yield.refEff));

    //Calculate the signal yield from reference yield, efficiencies and the Jpsi->mumu BR
    yield.sigYield = BR_sig / BR_ref * yield.refYield * yield.sigEff/yield.refEff;

    //Get the background yield from side band
    TChain * tree   = new TChain("DecayTree");

    std::vector<string> years;
    if (Run == 0) years.push_back(year);
    else years = yearsData(Run);

    for (auto& yr : years){
        string BDToutputPath = GetBDToutputFile(yr,getRunID(yr),false,false,false,KshortDecaysInVelo,UseLowQ2Range,false);
        tree->Add(BDToutputPath.c_str());
        coutDebug("Opening " + BDToutputPath);
    }

    string sVariable = (UseDTF ? "B_plus_M_DTF" : "B_plus_M");
    string cutMass = Form("%s > %f && %s < %f",
                           sVariable.c_str(), getBplusMeanFromResult(fitFile) -  effSigma, sVariable.c_str(),
                           getBplusMeanFromResult(fitFile) + effSigma);
    coutDebug("B mass cut: " + cutMass);

    //cut away Jpsi
    string cut = getFinalCut(false, false, false, "", gammaTMdefault, true, false, false, false, false, TMefficiencyClass(), -1, TMVAcut, RemoveMultiple);
    coutDebug("Using cut " + cut);

    //Get background estimation from sideband
    string bkgFitPath = GetBDTScanBackgroundFitFile(year,Run,KshortDecaysInVelo,UseLowQ2Range,TMVAcut);
    yield.bkgYield = GetBackgroundFromSidebandFit(tree, fixedMassRegion ? effSigma : B_plus_M_signal_window, cut, bkgFitPath);
    cut = cut + " && (" + cutMass + ")";
    coutDebug("Using cut " + cut);

    //data: create histograms from tree
    tree->Draw(Form("%s>>B_plus_M_plot", sVariable.c_str()), cut.c_str());
    TH1 * histo = ((TH1 *) gPad->GetPrimitive("B_plus_M_plot"));
    Int_t EventsAfterPreSelection = histo->GetEntries();

    yield.allSigEvts = EventsAfterPreSelection;

    //Print signal, background and all evts
    string trackType = (Kst2Kspiplus && SplitDDandLL ? (KshortDecaysInVelo ? "LL Tracks: " : "DD Tracks: ") : ":");
    coutDebug(TheDecay + "[SIGNAL]: "  + trackType + to_string(yield.sigYield));
    coutDebug(TheDecay + "[BCKGND]: "  + trackType + to_string(yield.bkgYield));
    coutDebug(TheDecay + "[ALLEVTS]: " + trackType + to_string(yield.allSigEvts));

    coutInfo("Expected yield is: " + to_string(yield.sigYield) + "\t+-" + to_string(yield.sigYieldErr));

    return yield;
}

YieldInfo GetSigAndBkgEstimationFromData(string year = "2011", int Run = 0, int randomSubset = 0,
                                         bool KshortDecaysInVelo = false, bool UseLowQ2Range = false,
                                         Double_t TMVAcut = -1.0, bool scan = false){
    gStyle->SetOptStat(0);
    LHCbStyle();

    bool fixedMassRegion = !Kst2Kspiplus; //Fix it in pi0 channel //TODO
    string magnet = "both"; //prepare in case of polarity studies

    //initialize the struct
    YieldInfo yield = YieldInfo();
    yield.year = year;

    //Set signal and background function separately for each decay channel
    string SignalFunction = "OneCB";
    string BackGroundFunction = GetBackgroundFunction(KshortDecaysInVelo);

    if (!scan) { //Do the fit first
        yield.refYield = massFit(year,magnet,Run,
                                 false, true, false, false,
                                 true, false,
                                 FixShape, SignalFunction, BackGroundFunction, FixShape,
                                 KshortDecaysInVelo, UseLowQ2Range,
                                 TMVAcut, randomSubset,
                                 fixedMassRegion, false,
                                 false, true, false,
                                 "", -1,
                                 RemoveMultiple, true,
                                 false, "", false,
                                 "", gammaTMdefault, false);
        yield.sigYield = massFit(year,magnet,Run,
                                 false, true, false, false,
                                 false, true,
                                 FixShape, SignalFunction, BackGroundFunction, FixShape,
                                 KshortDecaysInVelo, UseLowQ2Range,
                                 TMVAcut, randomSubset,
                                 fixedMassRegion, false,
                                 false, true, false,
                                 "", -1,
                                 RemoveMultiple, true,
                                 false, "", false,
                                 "", gammaTMdefault, false);
    }    
    //Load fit results now:

    // Get reference yield
    TFile *fitFileRef  = new TFile(GetMassFitFile(year, magnet, Run,
                                                  false, true, false, false,
                                                  true, false,
                                                  FixShape, SignalFunction, BackGroundFunction, FixShape,
                                                  KshortDecaysInVelo, UseLowQ2Range,
                                                  TMVAcut, randomSubset,
                                                  fixedMassRegion, false,
                                                  "", -1,
                                                  RemoveMultiple, false, //Data cannot be weighted technically, so the name is without "weighted"
                                                  false, "", false,
                                                  "", gammaTMdefault, false).c_str(),"OPEN");

    coutDebug("Opening " + string(fitFileRef->GetPath()));
    yield.refYield    = getSigYield(fitFileRef);
    yield.refYieldErr = getSigYieldErr(fitFileRef);
    fitFileRef->Close();

    // Get signal yield
    TFile *fitFile  = new TFile(GetMassFitFile(year, magnet, Run,
                                               false, true, false, false,
                                               false, true,
                                               FixShape, SignalFunction, BackGroundFunction, FixShape,
                                               KshortDecaysInVelo, UseLowQ2Range,
                                               TMVAcut, randomSubset,
                                               fixedMassRegion, false,
                                               "", -1,
                                               RemoveMultiple, false, //Data cannot be weighted technically, so the name is without "weighted"
                                               false, "", false,
                                               "", gammaTMdefault, false).c_str(),"OPEN");

    coutDebug("Opening " + string(fitFile->GetPath()));
    yield.sigYield    = getSigYield(fitFile);
    yield.sigYieldErr = getSigYieldErr(fitFile);
    yield.bkgYield    = getBkgYield(fitFile);
    yield.bkgYieldErr = getBkgYieldErr(fitFile);
    double mean_sig = getBplusMeanFromResult(fitFile);
    double effSigma = fixedMassRegion ? B_plus_M_signal_window : getEffSigma(fitFile);
    fitFile->Close();

    //Get the number of all events in the signal region
    TChain * tree   = new TChain("DecayTree");
    std::vector<string> years;
    if (Run == 0) years.push_back(year);
    else years = yearsData(Run);

    for (auto& yr : years){
        string BDToutputPath = GetBDToutputFile(yr,getRunID(yr),false,false,false,KshortDecaysInVelo,UseLowQ2Range,false);
        tree->Add(BDToutputPath.c_str());
        coutDebug("Opening " + BDToutputPath);
    }

    coutDebug("eff sigma: " + to_string(effSigma));
    string sVariable = (UseDTF ? "B_plus_M_DTF" : "B_plus_M");
    string cutMass = Form("%s > %f && %s < %f",
                           sVariable.c_str(), mean_sig -  effSigma, sVariable.c_str(),
                          mean_sig + effSigma);
    coutDebug("B mass cut: " + cutMass);

    //cut away resonances
    string cut = getFinalCut(false, false, false, "", gammaTMdefault, true, false, false, false, false, TMefficiencyClass(), -1, TMVAcut, RemoveMultiple);
    coutDebug("Using cut " + cut);

    //data: create histograms from tree
    tree->Draw(Form("%s>>B_plus_M_plot", sVariable.c_str()), cut.c_str());
    TH1 * histo = ((TH1 *) gPad->GetPrimitive("B_plus_M_plot"));
    Int_t EventsAfterPreSelection = histo->GetEntries();

    yield.allSigEvts = EventsAfterPreSelection;

    //Print signal, background and all evts
    string trackType = (Kst2Kspiplus && SplitDDandLL ? (KshortDecaysInVelo ? "LL Tracks: " : "DD Tracks: ") : ":");
    coutDebug(TheDecay + "[SIGNAL]: "  + trackType + to_string(yield.sigYield));
    coutDebug(TheDecay + "[BCKGND]: "  + trackType + to_string(yield.bkgYield));
    coutDebug(TheDecay + "[ALLEVTS]: " + trackType + to_string(yield.allSigEvts));

    coutInfo("Expected yield is: " + to_string(yield.sigYield) + "\t+-" + to_string(yield.sigYieldErr));
    coutDebug("Signal yield error: " + to_string(yield.sigYieldErr));
    coutDebug("Background yield error: " + to_string(yield.bkgYieldErr));
    coutDebug("Reference channel  yield error: " + to_string(yield.refYieldErr));

    if (yield.refYield == 0) coutERROR("Reference channel yield is zero!");
    if (yield.sigYield == 0) coutERROR("Signal channel yield is zero! Abort.");

    return yield;

}

int SaveTGraphs(string path, bool fineScan, string year, int Run, string basicPath,
                TGraphErrors *yieldGraph, TGraphErrors *bkgGraph, TGraphErrors *bkgGraphFromAllEvts, TGraphErrors *refYieldGraph, TGraphErrors *allEvtsInSig, TGraphErrors *significance, TGraphErrors *significanceFromAllEvts){

    coutInfo("Making pretty plots and saving them!");
    coutDebug("Writting into " + path);
    coutDebug("Using basicPath " + basicPath);

    TFile * TGraphOutput = new TFile(path.c_str(), "RECREATE");
    TGraphOutput->cd();

    designYieldGraph(yieldGraph,              Run, year, "sigYield",                 basicPath, TGraphOutput, false);
    designYieldGraph(bkgGraph,                Run, year, "bkgYield",                 basicPath, TGraphOutput, false);
    designYieldGraph(bkgGraphFromAllEvts,     Run, year, "bkgYield_fromAllEvts",     basicPath, TGraphOutput, fineScan);
    designYieldGraph(refYieldGraph,           Run, year, "refYield",                 basicPath, TGraphOutput, fineScan);
    designYieldGraph(significance,            Run, year, "significance",             basicPath, TGraphOutput, false);
    designYieldGraph(significanceFromAllEvts, Run, year, "significance_fromAllEvts", basicPath, TGraphOutput, fineScan);
    designYieldGraph(allEvtsInSig,            Run, year, "EventsInSigRegion",        basicPath, TGraphOutput, false);

    TGraphOutput->Close();
    return 1;


}

int ScanSignalAndBckgndEstimation(string year, int Run, Double_t BDTstep, bool KshortDecaysInVelo, bool UseLowQ2Range, bool scan, bool fineScan){
    //If Run==0, do the scan for year, if Run == 1/2, do the scan per run

    if (verboseLevel > 1)RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);

    bool section = true; //In case of fine scan, section the scan in more fine way
    int n_divisions = 3;

    //Get low TMVA response lower boundary
    Double_t lowBDTcut = (TMVAmethod == "MLP") ? 0.0 : -1.0;
    if (fineScan){
        coutInfo("Making fine BDT cut study.");
        lowBDTcut = 0.9;
        BDTstep =   0.001;
    }

    //prepare in case of polarity studies
    string magnet = "both";

    TGraphErrors *yieldGraph = new TGraphErrors();
    TGraphErrors *bkgGraph = new TGraphErrors();
    TGraphErrors *bkgGraphFromAllEvts = new TGraphErrors();
    TGraphErrors *refYieldGraph = new TGraphErrors();
    TGraphErrors *significance = new TGraphErrors();
    TGraphErrors *significanceFromAllEvts = new TGraphErrors();
    YieldInfo yield = YieldInfo(); //put zeroes everywhere

    //loop over BDT cuts
    coutInfo("Loop over BDT cuts with a step of " + to_string(BDTstep));
    coutDebug("Starting the loop at " + to_string(lowBDTcut));
    double tmpBDTcut = 0.0;

    for(double fBDTcut = lowBDTcut; fBDTcut < (section ? 0.999 : 1.0); fBDTcut += section ? ((1.0-fBDTcut)/double(n_divisions)) : BDTstep){
        for (int k = 0; k < n_divisions; k++){
            tmpBDTcut = fBDTcut + k*(1.0-fBDTcut)/double(n_divisions);
            yield = GetSigAndBkgEstimation(year, Run, KshortDecaysInVelo, UseLowQ2Range, tmpBDTcut, scan);
            yieldGraph              ->SetPoint(yieldGraph->GetN(),          tmpBDTcut, yield.sigYield);
            bkgGraphFromAllEvts     ->SetPoint(bkgGraphFromAllEvts->GetN(), tmpBDTcut, yield.allSigEvts-yield.sigYield);
            bkgGraph                ->SetPoint(bkgGraph->GetN(),            tmpBDTcut, yield.bkgYield);
            refYieldGraph           ->SetPoint(yieldGraph->GetN(),          tmpBDTcut, yield.refYield);
            significance            ->SetPoint(significance->GetN(),        tmpBDTcut, yield.sigYield / (TMath::Sqrt(yield.sigYield + yield.bkgYield)));
            significanceFromAllEvts ->SetPoint(significance->GetN(),        tmpBDTcut, yield.sigYield / (TMath::Sqrt(yield.allSigEvts)));
            coutDebug("---" + string(TMVAmethod) + " cut " + to_string(tmpBDTcut) + " yield: " + to_string(yield.sigYield) + ", background: " + to_string( yield.allSigEvts-yield.sigYield ));
        }
    }

    //Write to file and close
    string path = GetBDTScanFile(year, magnet, Run, KshortDecaysInVelo, UseLowQ2Range, fineScan);
    SaveTGraphs(path,fineScan,year,Run,path,yieldGraph,bkgGraph,bkgGraphFromAllEvts,refYieldGraph,NULL,significance,significanceFromAllEvts);

    coutInfo("Signal scan done for year " + year + "!");
    return 1;

}

int ScanSignalAndBckgndEstimationPerYear(string year = "2011", Double_t BDTstep = 0.01, bool KshortDecaysInVelo = true, bool UseLowQ2Range = false, bool scan = false, bool fineScan = false){
    return ScanSignalAndBckgndEstimation(year = "2011", 0, BDTstep, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
}
int ScanSignalAndBckgndEstimationAllYears(Double_t BDTstep = 0.01, bool KshortDecaysInVelo = true, bool UseLowQ2Range = false, bool scan = false, bool fineScan = true){
    for (auto yr: yearsData(12)) ScanSignalAndBckgndEstimationPerYear(yr, BDTstep, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
    return 1;
}

int ScanSignalAndBckgndEstimationSimple(string year = "2011",int Run = 1, Double_t BDTstep = 0.01, int randomSubset = 0,
                                        bool KshortDecaysInVelo = true, bool UseLowQ2Range = false, bool scan = false, bool fineScan = true){
    //Get low TMVA response lower boundary

    double lowBDTcut = (TMVAmethod == "MLP") ? 0.0 : -1.0;
    if (fineScan){
        coutInfo("Making fine BDT cut study.");
        lowBDTcut = 0.9;
        BDTstep = 0.002;
    }
    bool halve = true;


    //prepare in case of polarity studies
    string magnet = "both";

    TGraphErrors *yieldGraph = new TGraphErrors();
    TGraphErrors *bkgGraph = new TGraphErrors();
    TGraphErrors *bkgGraphFromAllEvts = new TGraphErrors();
    TGraphErrors *refYieldGraph = new TGraphErrors();
    TGraphErrors *allEvtsInSig = new TGraphErrors();
    TGraphErrors *significance = new TGraphErrors();
    TGraphErrors *significanceFromAllEvts = new TGraphErrors();
    YieldInfo yield = YieldInfo(); //put zeroes everywhere

    //loop over BDT cuts
    coutInfo("Loop over BDT cuts with a step of " + to_string(BDTstep));
    coutDebug("Starting the loop at " + to_string(lowBDTcut));
    int n_divisions = halve ? 3 : 1;
    double tmpcut = 0.0;
    for(double fBDTcut = lowBDTcut; fBDTcut < (halve ? 0.9999 : 1); fBDTcut += halve ? ((1.0-fBDTcut)/double(n_divisions)) : BDTstep){
        for (int k = 0; k < n_divisions-1; k++){
            tmpcut = fBDTcut + k*(1.0-fBDTcut)/double(n_divisions);
            coutDebug("Cutting at: " + to_string(tmpcut));
            yield =GetSigAndBkgEstimationFromData(year,Run, randomSubset, KshortDecaysInVelo, UseLowQ2Range, tmpcut, scan);
            //Save the efficiency
            yieldGraph              ->SetPoint(yieldGraph->GetN(),          tmpcut, yield.sigYield);
            yieldGraph              ->SetPointError(yieldGraph->GetN()-1,   0,       yield.sigYieldErr);
            bkgGraph                ->SetPoint(bkgGraph->GetN(),            tmpcut, yield.bkgYield);
            bkgGraph                ->SetPointError(bkgGraph->GetN()-1,     0,       yield.bkgYieldErr);
            bkgGraphFromAllEvts     ->SetPoint(bkgGraph->GetN(),            tmpcut, yield.allSigEvts-yield.sigYield);
            refYieldGraph           ->SetPoint(yieldGraph->GetN(),          tmpcut, yield.refYield);
            refYieldGraph           ->SetPointError(yieldGraph->GetN()-1,   0,       yield.refYieldErr);
            allEvtsInSig            ->SetPoint(allEvtsInSig->GetN(),        tmpcut, yield.allSigEvts);
            significance            ->SetPoint(significance->GetN(),        tmpcut, yield.sigYield / (TMath::Sqrt(yield.sigYield+yield.bkgYield)));
            significanceFromAllEvts ->SetPoint(significance->GetN(),        tmpcut, yield.sigYield / (TMath::Sqrt(yield.allSigEvts)));
            coutDebug("---" + string(TMVAmethod) + " cut " + to_string(tmpcut) + " yield: " + to_string(yield.sigYield) + ", background: " + to_string( yield.allSigEvts-yield.sigYield ));
        }
    }

    //Write to file and close //TODO at some point
    string path = GetBDTScanFile(year, magnet, Run, KshortDecaysInVelo, UseLowQ2Range, fineScan);
    replace(path, ".root", "fromData.root");
    if (randomSubset == -1) replace(path,".root","_subset1.root");
    if (randomSubset ==  1) replace(path,".root","_subset2.root");
    SaveTGraphs(path,fineScan,year,0,path,yieldGraph,bkgGraph,bkgGraphFromAllEvts,refYieldGraph,NULL,significance,significanceFromAllEvts);

    coutInfo("Signal scan done for Run " + to_string(Run) + "!");
    return 1;
}

int ScanSignalAndBckgndEstimationSimplePerYear(string year = "2011", Double_t BDTstep = 0.01, int randomSubset = 0, bool KshortDecaysInVelo = true, bool UseLowQ2Range = false, bool scan = false, bool fineScan = true){
    return ScanSignalAndBckgndEstimationSimple(year = "2011", 0, BDTstep, randomSubset, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
}

int ScanSignalAndBckgndEstimationSimpleAllYears(Double_t BDTstep = 0.01, int randomSubset = 0, bool KshortDecaysInVelo = true, bool UseLowQ2Range = false, bool scan = false, bool fineScan = true){
    for (auto yr: yearsData(12)) ScanSignalAndBckgndEstimationSimplePerYear(yr, BDTstep, randomSubset, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
    return 1;
}

double getMaxBDTresponse(string year = "2011", int Run = 0,  bool fineScan = false, bool directScan = false, int randomSubset = 0, bool KshortDecaysInVelo = true, bool UseLowQ2Range  = false){

    //prepare in case of polarity studies
    string magnet = "both";

    if (Run == 0) coutInfo("Getting BDTresponse with highest significance for YEAR " + year + ". .  .");
    else          coutInfo("Getting BDTresponse with highest significance for RUN "+to_string(Run)+". .  .");

    string path = GetBDTScanFile(year, magnet, Run, KshortDecaysInVelo, UseLowQ2Range, fineScan);
    if (directScan) replace(path, ".root", "fromData.root");
    if (randomSubset == -1) replace(path,".root","_subset1.root");
    if (randomSubset ==  1) replace(path,".root","_subset2.root");

    TGraphErrors *significance = new TGraphErrors();
    TFile *scanFile = new TFile(path.c_str(),"OPEN");
    string graphName = "significance_fromAllEvts";
    significance = (TGraphErrors*) scanFile->Get(graphName.c_str());
    if(significance == NULL){
        coutERROR("TGraphError was not found in the file!");
        coutERROR("File:\t" + path);
        coutERROR("TGraph:\t" + graphName);
        return 0;
    }

    //Get the best BDTcut
    double LargestFoM = -100.;
    double bestCut = 0.0;
    for(int p = 0; p < significance->GetN(); p++){
        if(significance->GetY()[p] > LargestFoM){
            LargestFoM  = significance->GetY()[p];
            bestCut     = significance->GetX()[p];
        }
    }
    return bestCut;
}

int optimizeBDTCut(string year = "2011", int Run = 0, bool KshortDecaysInVelo = true, Double_t BDTstep = 0.01, bool UseRandomSubset = false,
                       bool UseLowQ2Range  = false, bool scan = false, bool fineScan = false, bool directScan = false){

    //Get low TMVA response lower boundary
    double lowBDTcut = (TMVAmethod == "MLP") ? 0.0 : -1.0;
    if (fineScan){
        coutInfo("Making fine BDT cut study.");
        lowBDTcut = 0.9;
        BDTstep = 0.002;
    }

    double BestCut = 0;

    if (directScan){
        if (UseRandomSubset){
            ScanSignalAndBckgndEstimationSimple(year, Run, BDTstep, -1, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
            ScanSignalAndBckgndEstimationSimple(year, Run, BDTstep, +1, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
            double BestCut1 = getMaxBDTresponse(year, Run, fineScan, true, -1, KshortDecaysInVelo, UseLowQ2Range);
            double BestCut2 = getMaxBDTresponse(year, Run, fineScan, true, +1, KshortDecaysInVelo, UseLowQ2Range);
            BestCut = roundf(BestCut1 * 100) / 100 +roundf(BestCut2 * 100); // round to 2 decimal points, return e.g. 0.96 + 94 (easy to untangle)
        }
        else{
            ScanSignalAndBckgndEstimationSimple(year, Run, BDTstep, 0, KshortDecaysInVelo, UseLowQ2Range, scan, fineScan);
            BestCut = getMaxBDTresponse(year, Run, fineScan, true, 0, KshortDecaysInVelo,  UseLowQ2Range);
        }

        //TODO
    }
    else{ //TODO: add subset fits
        if (Run == 0){
            ScanSignalAndBckgndEstimationPerYear(year,BDTstep,KshortDecaysInVelo,UseLowQ2Range,scan,fineScan);
        }
        else{
            ScanSignalAndBckgndEstimation(year,Run,BDTstep,KshortDecaysInVelo,UseLowQ2Range,scan,fineScan);
        }
        //TODO make a pretty plot

        BestCut = getMaxBDTresponse(year, Run, fineScan, false, 0, KshortDecaysInVelo, UseLowQ2Range);

    }

    return BestCut;

}