kopecna/EWP-BplusToKstMuMu-AngAna
kopecna/EWP-BplusToKstMuMu-AngAna
2
1
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
e154e70c18
EWP-BplusToKstMuMu-AngAna/Code/Selection/MassFit.cpp

1054 lines
46 KiB
C++
Raw Normal View History

First commit, not including the most of the Codes/ folder
2022-02-10 15:49:54 +01:00
//Functions used for the B plus mass fit
//Renata Kopecna
#include "GlobalFunctions.hh"
#include "MassFit.hpp"
#include "BmassShape/BackgroundPdf.hpp"
#include "BmassShape/SignalPdf.hpp"
#include "Paths.hpp"
#include "MVAclass.hpp"
#include "RooFit/RooDoubleCB/RooDoubleCB.h"
#include "Design.hpp"
#include "Utils.hpp"
using namespace std;
using namespace RooFit;
using namespace RooStats;
bool useExtraVarBool(string extraVar){ //Shoudl be easy to add others
vector<string> allVars {"q2", "q2_binned", "q2_binned_fit", "thetak", "thetal", "phi", "pi_zero_ETA", "pi_zero_ETA-K_plus_ETA", "pi_zero_P", "pi_zero_P_DTF", "pi_zero_PT", "pi_zero_PT_DTF","K_plus_P","K_plus_PT"};
for (auto &var: allVars){
if (extraVar == var) return true;
if (extraVar == var+"_equal") return true;
}
return false;
}
string GetsWeightPlots(string year, bool UseOnlyJpsiEvents, bool UseOnlyMuMuEvents, bool KshortDecaysInVelo, bool GetShapeFromMC, string SignalType, string BkgType, bool ConstrainParameters){
string path = GetMassFitFile(year,false,UseOnlyJpsiEvents,UseOnlyMuMuEvents, KshortDecaysInVelo, GetShapeFromMC, SignalType,BkgType, ConstrainParameters,true,false);
replace(path,"BplusMassModel","Sweight");
return path;
}
double massFit(string year, string magnet, int Run,
bool MC, bool Preselected, bool TM, bool PHSP, //input/output file selection
bool UseOnlyJpsiEvents, bool UseOnlyMuMuEvents, //signal/reference
bool GetShapeFromMC, string SigType, string BkgType, bool ConstrainParameters, //shape
bool KshortDecaysInVelo, bool UseLowQ2Range, //LL/DD? q2range?
Double_t TMVAcut, int randomSubset, //TMVA options
bool fixedMassRegion, bool yieldOverFullRange, //yield calculation region
bool sWeight, //sWeight data?
bool loopFit, bool IsEfficiency, //additional options
string sExtraVar, int nExtraBin, //fit in bins of extra variable
bool removeMultiple, //Remove multiple candidates?
bool weighted, bool weightedFromPi0, string whichWeight, //use weight in the fit?
bool nonTM, string customTMbranch, bool gammaTM, //TM options
bool InclusiveSample){ //Fit inclusive sample?
//bool notTM: if true and MC, fits the nonTM part of MC with a CB
// To fit the CB as a background either in MC or Data, add "OneCB" to the BkgType string
//randomSubset \in {-1,0,1}
//loopFit means that MC fit is not performed, but parameters are fixed to previous fit, useful since one doesn't have to fit MC all the time in a loop
//IsEfficiency: Just stored in a different folder
//TODO: possibly add splittig in low/high Q2
//Make ROOT quiet
gStyle -> SetOptStat(0);
LHCbStyle();
gROOT->SetBatch(kTRUE);
//Load all libraries for ROOT
gSystem->Load("./BmassShape/SignalType_cpp.so");
gSystem->Load("./BmassShape/BackgroundType_cpp.so");
gSystem->Load("./BmassShape/ParamValues_cpp.so");
gSystem->Load("./BmassShape/SignalPdf_cpp.so");
gSystem->Load("./BmassShape/BackgroundPdf_cpp.so");
int printLevel = 1;
if (verboseLevel > 1){
printLevel = -1;
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
}
//sanity checks
magnet = correct_magnet_string(magnet);
checkMC(MC,false,PHSP,false);
if (!checkTM(MC,TM,nonTM,Preselected)) return 0;
checkKshort(KshortDecaysInVelo);
if (!SplitInQ2 && UseLowQ2Range) {
coutWarning("Cannot use low Q2 region if the data is not split in Q2! Setting UseLowQ2Range to false!");
UseLowQ2Range = false;
}
if (!checkQ2Range(UseOnlyJpsiEvents,UseOnlyMuMuEvents)) return 0;
if (TMVAcut > -1.0 && !Preselected){ //if TMVAcut == -1, no TMVA cut is applied
coutWarning("TMVA cut can be only aplied on preselected data! Setting Preselected to true!");
Preselected = true; //Cannot do BDT cut on stripped data
}
if(TM && MC && (GetShapeFromMC || ConstrainParameters)){
coutWarning("Cannot fix MC fit to an MC fit! Setting GetShapeFromMC and ConstrainParameters to false.");
GetShapeFromMC = false; //never run MC twice!
ConstrainParameters = false; // never use constrains on MC samples
}
if (sWeight && MC){
coutWarning("MC cannot be sWeighted! Setting sWeight to false!");
sWeight = false;
}
if (sWeight && UseOnlyMuMuEvents){
coutWarning("Data cannot be sWeighted only in the MuMu region! Setting UseOnlyMuMuEvents to false!");
UseOnlyMuMuEvents = false;
}
if (customTMbranch =="") customTMbranch = TMbranch;
if (customTMbranch =="TMedBKGCAT" && gammaTM == false) gammaTM = true;
bool useExtraVar = useExtraVarBool(sExtraVar); //perform once so it doesn't have to run through all the ifs every time
TMefficiencyClass extraVar = TMefficiencyClass(sExtraVar);
//To make life easier, make a bool for using RefChannel MC
bool UseRefMC = (MC && !PHSP && UseOnlyJpsiEvents);
//Make a years vecotr based on Run and selected proper RunID based on a year
std::vector<string> years;
int RunID = 0;
if (Run == 0){
RunID = getRunID(year);
if (MC && UseRefMC){
//can't use yearVector since we can replace 2017 and 2018 MC by 2016
if (!checkRefYear(year)) years.push_back("2016");
else years.push_back(year);
}
else years.push_back(checkIf2015MC(year,MC,UseRefMC,PHSP));
}
else{
RunID = Run;
if (InclusiveSample) years = yearsBkgMC(UseOnlyJpsiEvents,false, false,true, Run);
else if (MC) years = yearsMC(UseRefMC,PHSP,Run);
else years = yearsData(Run);
}
coutInfo("Using year/s ");
cout << "[INFO]\t"; printVector(years); cout << endl;
coutInfo(" with RunID "+to_string(RunID));
///signal shape
SignalPdf SigShape;
SigShape.setSignalType(SigType,GetShapeFromMC,ConstrainParameters);
///background shape
BackgroundPdf BkgShape;
BkgShape.setBackgroundType(BkgType,TM,GetShapeFromMC,ConstrainParameters);
///----------------------------------------------------
///Get Fit parameters from fit to TMed MC if requested:
///----------------------------------------------------
//Turn off MC fit in case of fitting in a loop, except for when TMVAcut == lowBDTcut
//When performing scans, always fix the parameters to the MC value
//however, the parameters are not calculated each iteration of the scan, because the MC shape does not change
Float_t lowBDTcutCheck = (TMVAmethod == "MLP") ? 0.0 : -1.0; //check if at lower limit
bool runMCforTMVAcut = (TMVAcut == lowBDTcutCheck || TMVAcut == 0.9); //if TMVA at lower limit/0.9 for fine scans, do the MC fit again
coutDebug("Run MC for TMVA cut? " +to_string(runMCforTMVAcut));
coutDebug("TMVA cut: " +to_string(TMVAcut));
bool MC_UseOnlyJpsiEvents, MC_UseOnlyMuMuEvents;
if(GetShapeFromMC && (runMCforTMVAcut || !loopFit)){
//In case of CB in the background model, fit the nonTMed MC by oneCB
if (sExtraVar == "q2_binned"){
MC_UseOnlyJpsiEvents = false;
MC_UseOnlyMuMuEvents = true;
}
else{
MC_UseOnlyJpsiEvents = UseOnlyJpsiEvents;
MC_UseOnlyMuMuEvents = UseOnlyMuMuEvents;
}
if (BkgShape.bkgOneCB){
std::cout << std::endl;
coutInfo("||==============================================||");
coutInfo("|| Fitting the MC to obtain background shape! ||");
coutInfo("||==============================================||");
std::cout << std::endl;
if(massFit(year,magnet,Run,
true,true,false,PHSP,
MC_UseOnlyJpsiEvents, MC_UseOnlyMuMuEvents,
false,"NoSig","OneCB",false,
KshortDecaysInVelo,UseLowQ2Range,
TMVAcut, 0, fixedMassRegion, yieldOverFullRange,
false, false, IsEfficiency, sExtraVar, nExtraBin,
removeMultiple, weighted, weightedFromPi0, whichWeight,
true, customTMbranch, gammaTMdefault, InclusiveSample) == 0){
coutERROR("Fit to MC data did not work out!");
return 0;
}
}
//Now fit the TMed MC
if (false){ //TODO: possibly switch to true
std::cout << std::endl;
coutInfo("||==================================================||");
coutInfo("|| Fitting the MVAed data to obtain signal shape! ||");
coutInfo("||==================================================||");
std::cout << std::endl;
if(massFit(year,magnet,Run,
false, true, false,PHSP,
true, false,
false, SigType,BkgType, false,
KshortDecaysInVelo, UseLowQ2Range,
getTMVAcut(RunID), 0, fixedMassRegion, yieldOverFullRange,
false, false, IsEfficiency, sExtraVar, nExtraBin,
true, false, false, "",
false, "", gammaTMdefault, false) == 0){
coutERROR("Fit to MC data did not work out!");
return 0;
}
}
else{
std::cout << std::endl;
coutInfo("||===========================================||");
coutInfo("|| Fitting the MC to obtain signal shape! ||");
coutInfo("||===========================================||");
std::cout << std::endl;
if(massFit(year,magnet,Run,
true,true, true,PHSP,
MC_UseOnlyJpsiEvents, MC_UseOnlyMuMuEvents,
false,SigType,"NoBckGnd",false,
KshortDecaysInVelo, UseLowQ2Range,
-1.0, 0, fixedMassRegion, yieldOverFullRange,
false, false, IsEfficiency, sExtraVar, nExtraBin,
removeMultiple, weighted, weightedFromPi0, whichWeight,
false, customTMbranch, gammaTM, InclusiveSample) == 0){
coutERROR("Fit to MC data did not work out!");
return 0;
}
}
std::cout << std::endl;
coutInfo("||===========================================||");
coutInfo("|| Apply signal shape and fit data with it! ||");
coutInfo("||===========================================||");
std::cout << std::endl;
}
//For now only MC can be weighted
if (weighted && !MC){
coutWarning("For now only weighted MC can be done! Setting weight to false!");
weighted = false;
}
if (weighted) SumW2Error(kTRUE);
//print info about signal and background shape
coutInfo("Using a " + SigType + " for the Signal shape!");
coutInfo("Using a " + BkgType + " for the Background shape!");
///------------------------------
///Load file
///------------------------------
TChain* tree = new TChain(treeName(MC,Preselected).c_str());
coutInfo("Reading data from TTree... ");
if (TMVAcut != -1 || weightedFromPi0 || removeMultiple) tree = get_BDT_TChain(years, MC, MC && UseOnlyJpsiEvents, PHSP, KshortDecaysInVelo, weightedFromPi0);
else if (weighted) tree = get_weighted_TChain(years, MC, UseOnlyJpsiEvents, PHSP, KshortDecaysInVelo);
else tree = get_basic_TChain(magnet,years,Preselected,MC,UseOnlyJpsiEvents,PHSP,false, false, InclusiveSample);
if(tree->GetEntries() == 0){
coutERROR("No entries found in TTree. Exit program!");
return 0;
}
else coutDebug("DONE!");
///------------------------------
/// Set branches
///------------------------------
string B_mass_branch = UseDTF ? (Preselected ? "B_plus_M_DTF" : "B_plus_DTF_M") : ("B_plus_M");
coutDebug("Deactivating all not needed branches.. ");
tree->SetBranchStatus("*",0);
//activating all needed branches
tree->SetBranchStatus("Q2", 1);
tree->SetBranchStatus(B_mass_branch.c_str(),1);
string weightBranch = weighted ? getWeightName(customTMbranch,gammaTM) : "";
if (weighted) coutDebug("WeightBranch = " + weightBranch + ".");
if (TMVAcut > -1.0) tree->SetBranchStatus(TMVAmethod+"response",1);
if (TM || nonTM) tree->SetBranchStatus(customTMbranch.c_str(),1);
if (TM || nonTM) tree->SetBranchStatus(gammaTMbranch.c_str(),1);
if (TM || nonTM) tree->SetBranchStatus(customTMbranch.c_str(),1);
if (randomSubset!=0) tree->SetBranchStatus("RandomSubSet",1);
if (Kst2Kspiplus && SplitDDandLL) tree->SetBranchStatus("KshortDecayInVeLo",1);
if (weighted) tree->SetBranchStatus(weightBranch.c_str(), 1);
if (weightedFromPi0) tree->SetBranchStatus(Form("w_%s",whichWeight.c_str()),1);
if (useExtraVar){
for (auto var:extraVar.sBranchName)
tree->SetBranchStatus(var.c_str(),1); //add last branch (or the only one)
}
if (removeMultiple) tree->SetBranchStatus(getAloneBranch(MC,TM,customTMbranch,gammaTM).c_str(),1);
coutInfo("DONE!");
///------------------------------
///Fill all needed variables in RooDataSet
///------------------------------
coutInfo("Loading the dataset...");
RooRealVar *B_plus_M;
B_plus_M = new RooRealVar(B_mass_branch.c_str(),
(Kst2Kspiplus ? (UseDTF ? "m(K_{S}^{0}#pi^{+}#mu^{+}#mu^{-})" : "m(#pi^{+}#pi^{-}#pi^{+}#mu^{+}#mu^{-})") : (UseDTF ? "m(K^{+}#pi^{0}#mu^{+}#mu^{-})" : "m(K^{+}#gamma#gamma#mu^{+}#mu^{-})") ),
get_cut_B_plus_M_low(year),
cut_B_plus_M_high,
" MeV/c^{2} ");
//In case we want to change binning, this is really not ideal but oh well
string finalCut = getFinalCut(MC,TM, nonTM, customTMbranch, gammaTM, UseOnlyMuMuEvents, UseOnlyJpsiEvents, SplitInQ2, UseLowQ2Range, useExtraVar, extraVar, nExtraBin, TMVAcut, removeMultiple);
if (!Preselected) finalCut = "B_plus_DTF_M >0"; //Dummy
coutDebug("Using cuts " + finalCut);
RooRealVar Q2Bin ("Q2", " ", -0.1e6, 21.1e6," ");
RooRealVar TMVAresponse (TMVAmethod+"response", " ", -1.0, 1.0, " ");
RooRealVar TMvar (customTMbranch.c_str(), " ", -0.1, 1.1, " ");
RooRealVar TMgammaVar (gammaTMbranch.c_str(), " ", -0.1, 6.1, " ");
RooRealVar KshortDecayInVeLo("KshortDecayInVeLo", " ", KshortDecaysInVelo ? 0.9 : -0.1, KshortDecaysInVelo ? 1.1 : 0.1," ");
RooRealVar RandomSubSet ("RandomSubSet", " ", randomSubset - 0.1, randomSubset + 0.1," ");
const int nExtraVars = extraVar.sBranchName.size();
RooRealVar *RooExtraVars[nExtraVars];
for_indexed (auto var:extraVar.sBranchName){
RooExtraVars[i] = new RooRealVar(var.c_str(), var.c_str(), extraVar.vVarRange.at(i).at(0), extraVar.vVarRange.at(i).at(1)," ");
coutDebug("Var: " + string(var.c_str()));
coutDebug("Range: " + to_string(extraVar.vVarRange.at(i).at(0)) + " to " + to_string(extraVar.vVarRange.at(i).at(1)));
}
RooRealVar IsAlone (getAloneBranch(MC,TM,customTMbranch,gammaTM).c_str(), " ", -1.0, 2.0, " ");
RooRealVar *FinalWeight;
if (weighted && weightedFromPi0){ //TODO: the question is if this works, test it at some point!
FinalWeight = new RooRealVar("weight","weight",-10.0,10.0,"");
double w;
double w_pi0;
tree->SetBranchAddress(weightBranch.c_str(),&w);
tree->SetBranchAddress(Form("w_%s",whichWeight.c_str()),&w_pi0);
for (int e = 0; e<tree->GetEntries();e++){
tree->GetEntry(e);
FinalWeight->setVal(w*w_pi0);
}
}
else if (weighted){
FinalWeight = new RooRealVar(weightBranch.c_str(), weightBranch.c_str(), -10.0, 10.0," ");
}
else if (weightedFromPi0){
FinalWeight = new RooRealVar(Form("w_%s",whichWeight.c_str()), Form("w_%s",whichWeight.c_str()), 0.0, 10.0," ");
}
RooArgSet argSet = RooArgSet(*B_plus_M);
argSet.add(Q2Bin);
if (TM||nonTM) argSet.add(TMvar);
if (TM||nonTM) argSet.add(TMgammaVar);
if (useExtraVar)
for(int i = 0; i < nExtraVars;i++) argSet.add(*RooExtraVars[i]);
if (TMVAcut > -1.0) argSet.add(TMVAresponse);
if (randomSubset) argSet.add(RandomSubSet);
if (SplitDDandLL&&Kst2Kspiplus) argSet.add(KshortDecayInVeLo);
if (removeMultiple) argSet.add(IsAlone);
if (weighted||weightedFromPi0) argSet.add(*FinalWeight);
RooDataSet* data= NULL;
if (weighted||weightedFromPi0) data = new RooDataSet("data","data",argSet,Import(*tree),Cut(finalCut.c_str()), WeightVar(*FinalWeight));
else data = new RooDataSet("data","data",argSet,Import(*tree),Cut(finalCut.c_str()));
coutInfo(" DONE!");
if (verboseLevel < 3) data->Print();
coutDebug("Using weighted events? " + to_string (data->isWeighted()));
coutDebug("The dataset contains " + to_string(data->numEntries()) + " events!");
///------------------------------
///Build the B+ mass signal shape
///------------------------------
gROOT->ProcessLine(".L ./RooFit/RooDoubleCB/RooDoubleCB_cpp.so");
gROOT->ProcessLine(".L ./GlobalFunctions.hh"); //This is very weird, but it works ¯\_(ツ)_/¯
coutInfo("Setting signal fit parameters ... ");
//set mean of B_plus mass peak
RooRealVar* meanBplus = SigShape.getRooRealVar("mean",true);
///------------------------------
///Yields and Background + Signal Fit-Model
///------------------------------
//yields
coutDebug("The dataset has " + to_string(data->numEntries()) + " entries!");
RooRealVar N_Bplus("N_Bplus", "#B^{+}", data->numEntries()/2., 0., 2*data->numEntries());
RooRealVar N_comb("N_comb","N_comb", data->numEntries()/2., 0., 2*data->numEntries());
//Create the signal PDF
//Check if one selected any sigPdf and bkgPdf
if (SigShape.NoSig && BkgShape.NoBackground){
coutERROR("No fit model was selected. Abort.");
return 0;
}
coutDebug("Creating BplusMassModel...");
RooAddPdf *BplusMassModel = NULL;
if (!SigShape.NoSig) BplusMassModel = SigShape.getBplusMassModel(B_plus_M,meanBplus, UseOnlyJpsiEvents);
coutInfo("Signal PDF is created!");
if (!SigShape.NoSig){
coutInfo("Constraints PDF:");
SigShape.ConsPDF->Print();
}
coutDebug("Creating BplusBckGndModel...");
///Build the B+ mass background model
RooAddPdf * BplusBckGndModel = NULL;
if (!BkgShape.NoBackground) BplusBckGndModel = BkgShape.getBplusBkgModel(B_plus_M);
coutInfo("Background PDF is created!");
//add all required pdfs
RooAbsPdf* pdf = NULL;
if (BkgShape.NoBackground){
//Create the signal PDF
pdf = new RooAddPdf("pdf", "pdf", RooArgList(*BplusMassModel), RooArgList(N_Bplus));
}
else{
if (SigShape.NoSig) pdf = new RooAddPdf("pdf", "pdf", RooArgList(*BplusBckGndModel), RooArgList(N_comb));
else pdf = new RooAddPdf("pdf", "pdf", RooArgList(*BplusMassModel, *BplusBckGndModel), RooArgList(N_Bplus, N_comb));
}
if (verboseLevel<2) pdf->Print();
///Gaussian Constraint to MC values
if (!SigShape.AtLeastOneConstrainFound) coutDebug("No parameters in signal costrained.");
if (!BkgShape.AtLeastOneConstrainFound) coutDebug("No parameters in background costrained.");
///------------------------------
///Fit
///------------------------------
RooFitResult *result = new RooFitResult("FitResult","FitResult");
//Add constrained PDFs
if (SigShape.AtLeastOneConstrainFound || BkgShape.AtLeastOneConstrainFound){
RooArgSet *constrPdfs;
if (SigShape.AtLeastOneConstrainFound) constrPdfs = new RooArgSet(*SigShape.ConsPDF);
if (BkgShape.AtLeastOneConstrainFound) constrPdfs = new RooArgSet(*BkgShape.BkgConsPDF);
RooProdPdf *FitPDF_constraint = new RooProdPdf("FitPDF_constraint","FitPDF_constraint",RooArgSet(*pdf,*constrPdfs));
if (!IsEfficiency){
coutDebug("Constrained parameters are:" );
SigShape.ConsParameter->Print();
coutDebug("The following PDF is fitted to the data:" );
FitPDF_constraint->Print();
coutDebug("Start the fit with constrains! " );
}
result = FitPDF_constraint->fitTo(*data,Save(kTRUE),Extended(kTRUE),NumCPU(3),Constrain(*constrPdfs),PrintLevel(printLevel),SumW2Error(weighted ? kTRUE : kFALSE));
delete FitPDF_constraint;
}
else{
coutDebug("The following PDF is fitted to the data:" );
if (verboseLevel<2) pdf->Print();
coutDebug("Start the fit! " );
result = pdf->fitTo(*data,Save(kTRUE),Extended(kTRUE),NumCPU(3),PrintLevel(printLevel),SumW2Error(weighted ? kTRUE : kFALSE));
}
coutDebug("The fit result is --------------- " );
if (verboseLevel<2) result->Print();
//Saving the results in case of constrianing to MC
SigShape.setValuesAndErrors(true);
BkgShape.setValuesAndErrors(false);
///------------------------------
///calculate # (signal)background events in signal region
///------------------------------
coutDebug("The effective sigma is: " + to_string(SigShape.getEffSigmaFromResult()));
if (fixedMassRegion) B_plus_M->setRange("SigRange",meanBplus->getVal() - B_plus_M_signal_window, meanBplus->getVal() + B_plus_M_signal_window );
else B_plus_M->setRange("SigRange",meanBplus->getVal() - SignalRegionNsigma * SignalFitParameter.SigmaEff, meanBplus->getVal() + SignalRegionNsigma * SignalFitParameter.SigmaEff);
coutInfo("Determining the signal and background yield");
//Get signal yield
RooAbsReal* S_fr = NULL;
Double_t S = 0.0, S_error_test = 0.0, S_error = 0.0;
if (!SigShape.NoSig){
S_fr = BplusMassModel->createIntegral(*B_plus_M,NormSet(*B_plus_M),Range("SigRange"));
S = S_fr->getVal() * N_Bplus.getVal();
S_error_test = S_fr->getPropagatedError(*result) * N_Bplus.getVal();
S_error = S_fr->getVal() *N_Bplus.getError();
}
//Get background yield
RooAbsReal* B_fr = nullptr;
Double_t B = 0.0, B_error = 0.0, B_error_test = 0.0;
if(!BkgShape.NoBackground){
B_fr = BplusBckGndModel->createIntegral(*B_plus_M,NormSet(*B_plus_M),Range("SigRange"));
B = B_fr->getVal() * N_comb.getVal();
B_error_test = B_fr->getPropagatedError(*result) * N_comb.getVal();
B_error = B_fr->getVal() * N_comb.getError(); //B * (N_comb.getError()/N_comb.getVal());
}
coutInfo("S/sqrt(S+B)= " + to_string(S/sqrt(S+B)));
if(B != 0.)coutInfo("S/B= " + to_string(S/B));
coutInfo("S= " + to_string(S)+ "+-" + to_string(S_error));//+ "/+-" + to_string(S_error_test));
coutInfo("B= " + to_string(B)+ "+-" + to_string(B_error));//+ "/+-" + to_string(B_error_test));
cout << endl << endl;
///------------------------------
///Plot the models with indidivual components
///------------------------------
//configurables for plotting
Float_t pullHeight = 0.32;
Int_t Bins = 100;
if (UseOnlyMuMuEvents && Preselected && !MC) Bins = useExtraVar ? 30 : 50;
//create two pads into one canvas
TCanvas* c1= new TCanvas("");
TPad *pad1 = new TPad("pad1", "plot",0.0,pullHeight,1.0,1.0,0);
TPad *pad2 = new TPad("pad2", "pull",0.0,0.0,1.0,pullHeight,0);
pad1->Draw();
pad2->Draw();
//modifiy pad for main plot
pad1->SetBorderSize (0);
pad1->SetBottomMargin(1e-6);
pad1->SetTopMargin(pad1->GetTopMargin() / pullHeight * ( 1 - pullHeight) );
pad1->cd();
RooPlot* frame_m = B_plus_M->frame();
//some axis-modifications on the main plot
TGaxis::SetExponentOffset(1e+9,1e+9,"y");//offset = pad size * 1e+7
//make the frame pretty
designMassFitFrame(frame_m, pullHeight);
data->plotOn(frame_m,Name("data"),MarkerSize(0.5),Binning(Bins),DataError(RooAbsData::SumW2));
pdf->plotOn(frame_m,Name("pdf"),LineColor(kBlack),LineWidth(2));
if (!SigShape.NoSig) pdf->plotOn(frame_m,Components("BplusMassModel"),LineColor(kBlue),LineStyle(kDashed),LineWidth(1));
if (!BkgShape.NoBackground){pdf->plotOn(frame_m,Components("BplusBckGndModel"),LineColor(kRed),LineStyle(kDashed),LineWidth(1));
//in case you want to plot ExpGaus
//pdf->plotOn(frame_m,Components("bkg_exp1"),LineColor(kRed),LineStyle(kDashed),LineWidth(1));
//pdf->plotOn(frame_m,Components("ExpG"),LineColor(kRed),LineStyle(kDashed),LineWidth(1));
}
//Plot the pdf and the box with fit results into the frame
data->plotOn(frame_m,Name("data"),MarkerSize(0.5),Binning(Bins));
//pdf->paramOn(frame_m, Format("NEU",AutoPrecision(1)));
pdf->paramOn(frame_m, ShowConstants(kTRUE), Format("NEU",AutoPrecision(2)), Layout(0.65, 0.95, 0.88));
frame_m->Draw();
c1->Update();
// TODO: Add a box with constrained parameters posibly
//signal/background yield
Float_t posX = MC ? 0.40 : 0.41, posY = 0.84;
if(Kst2Kspiplus) posX = 0.34;
TLatex* fitresult = new TLatex();
fitresult->SetTextFont(132);
fitresult->SetTextColor(1);
fitresult->SetTextSize(0.04);
fitresult->SetTextAlign(13);
fitresult->SetNDC(1);
//TODO
if (!SigShape.NoSig) fitresult->DrawLatex(posX, posY - 0.08, Form("Signal: %.0f #pm %.0f", S, S_error)); //Signal defined by N sigma around peak mean value or in the mass window, depends on bool fixedMassRegion and B_plus_M_signal_window
if(!BkgShape.NoBackground) fitresult->DrawLatex(posX, posY - 0.13, Form("Background: %.0f #pm %.0f", B, B_error));
frame_m->addObject(fitresult) ;
frame_m->Draw("SAME");
if (MC) addLHCbtag(posX, 0.82, "simulation", 1./(1. - pullHeight));
else addLHCbtag(posX, 0.82, "data", 1./(1. - pullHeight));
//create pull histogram
pad2->Clear();
pad2->SetTopMargin(1e-6);
pad2->SetBottomMargin(pad2->GetBottomMargin() / pullHeight );
pad2->cd();
//Design pull histogram
RooPlot* pullFrame = B_plus_M->frame();
designPullFrame(pullFrame,frame_m,pullHeight);
//Add two 3-sigma lines:
TLine * lUp = threeSigmaLine(true);
TLine * lLow = threeSigmaLine(false);
pullFrame->addObject(lUp);
pullFrame->addObject(lLow);
//Draw the pull
pullFrame->Draw() ;
//choose path to save files //TODO: check
string filePath = "";
if (IsEfficiency) filePath = GetEfficiencyMassFitFile(year, magnet, Run,
Preselected, TM, PHSP, UseOnlyJpsiEvents, UseOnlyMuMuEvents,
GetShapeFromMC, SigType, BkgType, ConstrainParameters,
KshortDecaysInVelo, UseLowQ2Range,
TMVAcut, fixedMassRegion,
sExtraVar, nExtraBin,
removeMultiple, weighted, weightedFromPi0,
whichWeight, customTMbranch, gammaTM);
else filePath = GetMassFitFile(year,magnet,Run,MC,Preselected,TM,PHSP,UseOnlyJpsiEvents,UseOnlyMuMuEvents,
GetShapeFromMC, SigType, BkgType, ConstrainParameters,
KshortDecaysInVelo, UseLowQ2Range, TMVAcut, randomSubset,
fixedMassRegion, yieldOverFullRange, sExtraVar, nExtraBin,removeMultiple,
weighted, weightedFromPi0, whichWeight, nonTM, customTMbranch, gammaTM, InclusiveSample);
coutDebug("Saving into file "+ filePath);
TFile *fitFile = new TFile(filePath.c_str(),"RECREATE");
//save canvas to pdf and root format
coutInfo("Writting signal and background pdfs into the file " + filePath + ".");
fitFile->cd();
result->Write();
c1->Write();
//Save Signal yield, background yield and it's error
TVectorD yield(1);
yield[0] = S;
yield.Write("yield");
TVectorD yield_err(1);
yield_err[0] = S_error;
yield_err.Write("yield_err");
TVectorD B_vec(1);
B_vec[0] = B;
B_vec.Write("background");
TVectorD B_vec_err(1);
B_vec_err[0] = B_error;
B_vec_err.Write("background_err");
TVectorD sigma_eff(1);
sigma_eff[0] = SignalFitParameter.SigmaEff;
sigma_eff.Write("sigma_eff");
fitFile->Close();
coutDebug("File " + filePath + " closed." );
//Plot the fit
replace(filePath,".root",".eps");
coutInfo("Plotting into " + filePath + ".");
c1->Print(filePath.c_str());
if(sWeight){
if (BkgType == "NoBckGnd"){ //It crashes here if no background is selected
coutERROR("No background selected! Splot needs SingleExponential, DoubleExponential or ExpGauss.");
return 0;
}
meanBplus->setConstant();
SigShape.setAllRooVarsConstant();
BkgShape.setAllRooVarsConstant();
SPlot* sData = new SPlot("sData","An SPlot",*data,pdf,RooArgList(N_Bplus, N_comb));
///Plot the sWeight distributions as a function of mass
TCanvas* SwBplus = new TCanvas("Bplus sWeight","Bplus sWeight distribution");
TH2 * SwBplusHist;
SwBplusHist = (TH2*)data->createHistogram((UseDTF ? "B_plus_M_DTF,N_Bplus_sw" : "B_plus_M,N_Bplus_sw"));
SwBplusHist->GetYaxis()->SetTitle("Signal sWeights");
SwBplusHist->SetTitle("");
//SwBplus->Write("",TObject::kWriteDelete);
SwBplusHist->Draw();
addLHCbtag(0.65, 0.85, "data", 1);
filePath = GetsWeightPlots(year, UseOnlyJpsiEvents,UseOnlyMuMuEvents, KshortDecaysInVelo,GetShapeFromMC,SigType,BkgType,ConstrainParameters);
SwBplus->SaveAs(filePath.c_str());
replace(filePath,".root",".eps");
SwBplus->Print(filePath.c_str());
///Create output file
TFile* output = nullptr;
output = new TFile(GetBDTinputFile(year,false,UseOnlyJpsiEvents,PHSP,KshortDecaysInVelo).c_str() ,"RECREATE");
if(!output->IsOpen()){
coutERROR("Could not create output file. Abort!");
return 0;
}
output->cd();
tree->SetBranchStatus("*",1);
coutInfo("Copy Tree... " );
TTree * new_tree = nullptr;
string sVariable = (UseDTF ? "B_plus_M_DTF" : "B_plus_M");
string q2Cut = "";
if (UseOnlyJpsiEvents) //If for OnlyMuMu not needed so far //TODO
q2Cut = "&& " + getJpsicut();
if (SplitInQ2){
if (UseLowQ2Range)
q2Cut = "&& Q2 < 8.68e6";
else
q2Cut = "&& 10.09e6 < Q2";
}
//no use of saving only mumu data
if(Kst2Kspiplus && SplitDDandLL){
new_tree = tree->CopyTree(Form("%s >= %f && %s <= %f %s && KshortDecayInVeLo==%i", sVariable.c_str(), get_cut_B_plus_M_low(year), sVariable.c_str(), cut_B_plus_M_high, q2Cut.c_str(), KshortDecaysInVelo));
}
else{
new_tree = tree->CopyTree(Form("%s >= %f && %s <= %f %s", sVariable.c_str(), get_cut_B_plus_M_low(year), sVariable.c_str(), cut_B_plus_M_high, q2Cut.c_str()));
}
coutDebug("Finished!" );
double w;
TBranch* Bra_sw = new_tree->Branch("N_Bplus_sw", &w, "N_Bplus_sw/D");
///loop over events
int numEvents = new_tree->GetEntries();
coutDebug("Entries in TTree:\t"+ to_string(numEvents) );
coutDebug("Entries in RooDataSet:\t"+to_string(data->numEntries()));
if(numEvents != data->numEntries()){
coutERROR("Number of weights not equal to number of events");
return 0;
}
coutDebug("Loop over data sample " + year + TheDecay + " to save sWeights from mass fit results!" );
for(int i = 0; i < numEvents; i++){
if ((0ul == (i % 10000ul) || i + 1 == numEvents) && i != 0) coutInfo("Read event " + to_string(i) + "/" + to_string(numEvents));
tree->GetEntry(i);
w = sData->GetSWeight(i,"N_Bplus_sw");
Bra_sw->Fill();
}
coutInfo("Loop finished!!!");
output->cd();
new_tree->Write("",TObject::kWriteDelete);
delete new_tree;
output->Close();
delete SwBplus;
delete SwBplusHist;
delete sData;
}
if (S < 0.1) S = 0.01; //protect zeroes from happening
delete frame_m;
delete tree;
delete result;
coutInfo("Mass fit is done.");
return S;
}
int quickFit(string year, bool MC, bool sWeight, bool UseOnlyJpsiEvents, bool UseOnlyMuMuEvents, bool KshortDecaysInVelo, bool GetShapeFromMC, string SigType, string BkgType, bool ConstrainParameters){
return massFit(year,"both", 0,
MC, true, MC, false,
UseOnlyJpsiEvents, UseOnlyMuMuEvents,
GetShapeFromMC,SigType,BkgType,ConstrainParameters,
KshortDecaysInVelo,false,
-1.0, 0,
false, false,
sWeight,
false, false,
"",-1,
false,
false, false, "",
false, "", false,
false);
}
int quickTest(bool gammaTM = false, string customTMbranch = "TMed"){
setVerboseLevel(1);
return massFit("2015","both", 0,
false, true, true, false,
false, false,
true,"OneCB","SingleExponential",true,
false,false,
-1.0, 0,
false, false,
false,
false, false,
"",-1,false, true, false, "",false, customTMbranch, gammaTM,
false);
}
int testOneCbBackground(bool MC = false, string SigType = "OneCB", string BkgType = "SingleExponentialOneCB", string customTMbranch = "", bool gammaTM = false){
setVerboseLevel(1);
return massFit("2016","both",0,
false, true, true, false,
!MC, false,
true,SigType,BkgType,true,
false, false,
-1.0, 0,
true, false,
false, false, false,
"",-1, false, false, false, "", false, customTMbranch, gammaTM, false);
}
int efficiencyFit(string year = "2011", string magnet = "down", int Run = 0,
bool Preselected = true, bool TM = true, bool PHSP = false, //input/output file selection
bool UseOnlyJpsiEvents = false, bool UseOnlyMuMuEvents = false, //signal/reference
bool GetShapeFromMC = false, string SigType = "OneCB", string BkgType = "SingleExponential", bool ConstrainParameters = false, //shape
bool KshortDecaysInVelo = Kst2Kspiplus, //LL/DD?
Double_t TMVAcut = -1.0, //TMVA options
bool fixedMassRegion = !Kst2Kspiplus, //yield calculation region
bool UseLowQ2Range =false, //q2 ranges
string sExtraVar = "", int nExtraBin = -1,bool removeMultiple = false, bool weighted = false,
bool weightedFromPi0 = false, string whichWeight = "", //fit in bins of extra variable
string customTMbranch = "", bool gammaTM = false
){
return massFit(year,magnet,Run,
true, Preselected, TM, PHSP,
UseOnlyJpsiEvents, UseOnlyMuMuEvents,
GetShapeFromMC,SigType,BkgType,ConstrainParameters,
KshortDecaysInVelo, UseLowQ2Range,
TMVAcut, 0,
fixedMassRegion, false,
false, false, true,
sExtraVar, nExtraBin, removeMultiple,
weighted, weightedFromPi0,whichWeight,
false, customTMbranch, gammaTM, false
);
}
int basicYieldFit(string year, int Run, bool MC, bool PHSP, //input/output file selection
bool UseOnlyJpsiEvents, bool UseOnlyMuMuEvents, //signal/reference
bool GetShapeFromMC, string SigType, string BkgType, bool ConstrainParameters, //shape
bool KshortDecaysInVelo, bool UseLowQ2Range, //LL/DD? q2range?
Double_t TMVAcut, //TMVA options
bool fixedMassRegion, bool loopFit,//yield calculation region
bool removeMultiple //remove multiple candidates
){ //TODO: shape from MC if TMVA cut?
return massFit(year,"both",Run,
MC, true, MC, PHSP,//if MC then truthMatched
UseOnlyJpsiEvents, UseOnlyMuMuEvents,
GetShapeFromMC,SigType,BkgType,ConstrainParameters,
KshortDecaysInVelo, UseLowQ2Range,
TMVAcut, 0,
fixedMassRegion, false,
false, loopFit, false,
"", -1, removeMultiple,
true, false, "", //weighted=true only for MC, but also data then takes the shape from MC that is weighted
false, "", gammaTMdefault, false);
}
int basicYieldFitAllYears(//Fits all years separately
bool MC = true, bool PHSP = false, //input/output file selection
bool UseOnlyJpsiEvents = false, bool UseOnlyMuMuEvents = false, //signal/reference
bool GetShapeFromMC = false, string SigType = "OneCB", string BkgType = "SingleExponential", bool ConstrainParameters = false, //shape
bool KshortDecaysInVelo = Kst2Kspiplus, bool UseLowQ2Range =false, //LL/DD? q2range?
Double_t TMVAcut = -1.0, //TMVA options
bool fixedMassRegion = !Kst2Kspiplus, bool loopFit = false, //yield calculation region
bool removeMultiple = false //remove multiple candidates
){
for(auto &y : yearsVector(MC,UseOnlyJpsiEvents,false, 12)){
if (basicYieldFit(y,0, MC, PHSP, UseOnlyJpsiEvents,UseOnlyMuMuEvents,GetShapeFromMC, SigType, BkgType, ConstrainParameters,
KshortDecaysInVelo,UseLowQ2Range,TMVAcut, fixedMassRegion,loopFit, removeMultiple) ==0) return 0;
}
return 1;
}
int basicFitAllYearsAndRegions(bool MC = true, bool PHSP = false, //input/output file selection
bool GetShapeFromMC = false, string SigType = "OneCB", string BkgType = "SingleExponential", bool ConstrainParameters = false, //shape
bool KshortDecaysInVelo = Kst2Kspiplus, bool UseLowQ2Range =false, //LL/DD? q2range?
Double_t TMVAcut = -1.0, bool removeMultiple = false //TMVA options
){
vector <bool> UseOnlyJpsiEvents = {true,false};
vector <bool> UseOnlyMuMuEvents = {true,false};
vector <bool> fixedMassRegion = {true,false};
for(auto const&jpsi: UseOnlyJpsiEvents){
for(auto const&mumu: UseOnlyMuMuEvents){
for(auto const&fixedWindow: fixedMassRegion){
if (basicYieldFitAllYears(MC,PHSP, jpsi,mumu,GetShapeFromMC, SigType, BkgType, ConstrainParameters,
KshortDecaysInVelo,UseLowQ2Range,TMVAcut, fixedWindow,false,removeMultiple)==0) return 0;
}
}
}
return 1;
}
int getYieldBasicOptionsMC(bool TM, string year, int Run, bool ReferenceChannel, bool PHSP, bool removeMultiple){
string bkg = TM ? "NoBckGnd" : "SingleExponential";
bool useJpsiOnly = ReferenceChannel;
bool useMuMuOnly = !PHSP && !ReferenceChannel;
bool weighted = false;
return massFit(year,"both",Run,
true, true, TM, PHSP,
useJpsiOnly, useMuMuOnly,
false, "OneCB",bkg,false,
false, false,
-1.0, 0,
true, false,
false, false, false,
"", -1, removeMultiple,
weighted, false, "", //weighted=true only for MC, but also data then takes the shape from MC that is weighted
false, "", gammaTMdefault, false);
}
int getYieldBasicOptionsMCAllYears(bool TM, bool ReferenceChannel, bool PHSP, bool removeMultiple){
for (auto year: yearsMC(ReferenceChannel,PHSP,12)){
if (getYieldBasicOptionsMC(TM, year, 0, ReferenceChannel, PHSP, removeMultiple)==0) return 0;
}
if (getYieldBasicOptionsMC(TM, "2011", 1, ReferenceChannel, PHSP, removeMultiple) == 0) return 0;
if (getYieldBasicOptionsMC(TM, "2015", 2, ReferenceChannel, PHSP, removeMultiple) == 0) return 0;
return 1;
}
int getYieldAllYearsBasicOptions(bool TM, bool removeMultiple){
if (getYieldBasicOptionsMCAllYears(TM,false,false, removeMultiple)==0) return 0;
if (getYieldBasicOptionsMCAllYears(TM,true ,false, removeMultiple)==0) return 0;
if (getYieldBasicOptionsMCAllYears(TM,false,true, removeMultiple)==0) return 0;
//if (basicYieldFitAllYears(false, false, true, false, true, "OneCB", "SingleExponential", true,
// false, false, -1, true,false,true)==0) return 0;
//if (basicYieldFitAllYears(false, false, false, true, true, "OneCB", "SingleExponential", true,
// false, false, -1, true,false,true)==0) return 0;
return 1;
}
int basicYieldFitAllRuns( //Fits data per year
bool MC = true, bool PHSP = false, //input/output file selection
bool UseOnlyJpsiEvents = false, bool UseOnlyMuMuEvents = false, //signal/reference
bool GetShapeFromMC = false, string SigType = "OneCB", string BkgType = "SingleExponential", bool ConstrainParameters = false, //shape
bool KshortDecaysInVelo = Kst2Kspiplus, bool UseLowQ2Range =false, //LL/DD? q2range?
Double_t TMVAcut = -1.0, //TMVA options
bool fixedMassRegion = !Kst2Kspiplus, bool loopFit = false,//yield calculation region
bool removeMultiple = false //remove multiple candidates
){
if ( basicYieldFit("2011",1, MC, PHSP, UseOnlyJpsiEvents,UseOnlyMuMuEvents,GetShapeFromMC, SigType, BkgType, ConstrainParameters,
KshortDecaysInVelo,UseLowQ2Range,TMVAcut, fixedMassRegion,loopFit,removeMultiple)==0) return 0;
if ( basicYieldFit("2015",2, MC, PHSP,UseOnlyJpsiEvents,UseOnlyMuMuEvents,GetShapeFromMC, SigType, BkgType, ConstrainParameters,
KshortDecaysInVelo,UseLowQ2Range,TMVAcut, fixedMassRegion,loopFit,removeMultiple)==0) return 0;
return 1;
}
int massFitTestAll(bool GetShapeFromMC = true,
bool ConstrainParameters = false, int Run = 1){
bool KshortDecaysInVelo= false;
bool UseLowQ2Range = false;
bool TM = false;
bool MC = true;
bool UseOnlyJpsiEvents = false;
bool UseOnlyMuMuEvents = false;
bool UseFixedMassRegion = false;
bool yieldOverFullRange = false;
std::vector<string> years = yearsData(Run);
std::vector<string> signalShape = {//"SingleGaussian",
"DoubleGaussian",
//"CBLeft",
//"CBRight",
"CBDouble"
};
std::vector<string> backgroundShape = {//"SingleExponential",
"DoubleExponential",
"ExpGaus"
};
std::vector<string> polarity = { //"down",
//"up",
"both"
};
for (auto& y: years){
for (auto& sig: signalShape){
for (auto& bkg: backgroundShape){
for (auto& magnet: polarity){
if(massFit(y, magnet, 0, MC, true, TM, false,UseOnlyJpsiEvents, UseOnlyMuMuEvents, GetShapeFromMC, sig, bkg, ConstrainParameters,
KshortDecaysInVelo, UseLowQ2Range, -1.0, false, UseFixedMassRegion,yieldOverFullRange,false, false, false, "", -1,false, false, false, "",false,"",false,false) == 0)return 0;
}
}
}
}
return 1;
}
//check distributions in q2
int massFitTestQ2All(int Run, bool PHSP, bool UseOnlyMuMuEvents, bool UseOnlyJpsiEvents){
bool KshortDecaysInVelo= false;
bool UseLowQ2Range = false;
bool UseFixedMassRegion = false;
bool yieldOverFullRange = false;
bool removeMultiple = false;
bool weighted = PHSP;
string signalShape = "OneCB";
string backgroundShape = "ExpGaus";
for (int n = 4; n < 5; n++){
if(massFit("2011", "both", Run,
PHSP, true, PHSP, PHSP,
UseOnlyJpsiEvents, UseOnlyMuMuEvents,
true, signalShape, backgroundShape, true,
KshortDecaysInVelo, UseLowQ2Range,
-1.0, false, UseFixedMassRegion,yieldOverFullRange,
false, false, false, "q2_binned", n, removeMultiple, weighted,
false, "",false,
"",gammaTMdefault,false) == 0)return 0;
}
return 1;
}
//Print efficiencies and fit parameters
int PrintFitResults(RooFitResult* fitRes){
RooArgSet fitargs = fitRes->floatParsFinal();
if (fitargs.getSize() <=0){
coutERROR("RooArgSet size is zero/nonDef!");
return 0;
}
TIterator* iter(fitargs.createIterator());
cout << " ======= results ====== " << endl;
for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);
string name = rrv->GetName();
cout<< name << ": " << rrv->getVal() <<" +/- "<<rrv->getError() << endl;
}
return 1;
}
int fitJpsi(string year, int Run, bool MC, double TMVAcut, bool RemoveMultiple){
return massFit(year,"both", Run,
MC, true, MC, false,
true, false,
!MC, "OneCB", "SingleExponential", true,
false,false,
TMVAcut, 0,
false, false,
false,
false, false,
"",-1,
RemoveMultiple,
false, false, "",
false, "", gammaTMdefault,
false);
}
Reference in New Issue Copy Permalink