You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
636 lines
30 KiB
636 lines
30 KiB
//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;
|
|
|
|
}
|