diff --git a/Makefile b/Makefile index 1f3912e..f9b3c52 100644 --- a/Makefile +++ b/Makefile @@ -15,14 +15,13 @@ EXECUTABLE3 := simfit_genMC EXECUTABLE4 := simfit_genMC_multiFit EXECUTABLE5 := plotMultiResults EXECUTABLE6 := simfit_toy_fullAngular -EXECUTABLE7 := simfit_recoMC_fullAngularMass -EXECUTABLE8 := simfit_recoMC_fullMass -EXECUTABLE9 := simfit_recoMC_fullAngularMass_toybkg -EXECUTABLE10 := simfit_data_fullAngularMass -EXECUTABLE11 := simfit_data_fullAngularMass_Swave -EXECUTABLE12 := plot_simfit_data_fullAngularMass_Swave -EXECUTABLE13 := plot_simfit_recoMC_fullAngularMass_toybkg -EXECUTABLE14 := plot_simfit_recoMC_fullAngularMass + +file_with_mass_list = plot_ctk_slices_data plot_ctl_slices_data plot_phi_slices_data plot_mass_slices_data \ + plot_simfit_recoMC_fullAngularMass plot_simfit_recoMC_fullAngularMass_toybkg \ + plot_simfit_data_fullAngularMass_Swave simfit_data_fullAngularMass_Swave \ + simfit_data_fullAngularMass simfit_recoMC_fullAngularMass_toybkg \ + simfit_recoMC_fullMass simfit_recoMC_fullAngularMass \ + generateToy simfit_toy_fullAngularMass_Swave EXTRACLASS := RooDataHist.cxx CLASS0 := PdfRT @@ -82,31 +81,10 @@ $(EXECUTABLE5): $(EXECUTABLE5).cc $(EXECUTABLE6): $(EXECUTABLE6).cc $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) -$(EXECUTABLE7): $(EXECUTABLE7).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - -$(EXECUTABLE8): $(EXECUTABLE8).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - -$(EXECUTABLE9): $(EXECUTABLE9).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - -$(EXECUTABLE10): $(EXECUTABLE10).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) -$(EXECUTABLE11): $(EXECUTABLE11).cc +$(file_with_mass_list): %: %.cc $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - -$(EXECUTABLE12): $(EXECUTABLE12).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - -$(EXECUTABLE13): $(EXECUTABLE13).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - -$(EXECUTABLE14): $(EXECUTABLE14).cc - $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS) $(SOURCEDIR)/$(CLASS4).cc $(CLASSDICT2).cc $(ROOTLIBS) $(ROOTFLAGS) -I$(INCLUDEDIR) - #cleaning options .PHONY: clean clean: - rm -f $(EXECUTABLE0) $(EXECUTABLE1) $(EXECUTABLE7) + rm -f $(EXECUTABLE0) $(EXECUTABLE1) $(EXECUTABLE2) diff --git a/compareSimFitResultsResonantBins.cc b/compareSimFitResultsResonantBins.cc new file mode 100644 index 0000000..7e4897d --- /dev/null +++ b/compareSimFitResultsResonantBins.cc @@ -0,0 +1,319 @@ +using namespace RooFit; +using namespace std; + +/* +code to plot results from the Simultaneous fit of multiple years and compare them to +- GEN level results +- results of the fit to the single year datasets + (should be provided) +*/ + +int thebin = 6; + +static const int nBins = 8; // was 8 +float binBorders [nBins+1] = { 1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; + +static const int nPars = 8; +string ParName [nPars] = { "Fl", "P1", "P2", "P3", "P4p", "P5p", "P6p", "P8p" }; +string ParTitle [nPars] = { "F_{L}", "P_{1}", "P_{2}", "P_{3}", "P'_{4}", "P'_{5}", "P'_{6}", "P'_{8}" }; +float diffMax = 0.0999; + +int colorlist[5]; +TCanvas* c[nPars]; + +// for GEN +void fillVectors(RooFitResult* fitResult, double Res[nPars], double ErrH[nPars], double ErrL[nPars]){ + std::cout << "forGEN " << std::endl; + for (int ParIndx=0; ParIndxfloatParsFinal().find(ParName[ParIndx].c_str()); + Res[ParIndx] = Par->getValV(); + ErrH[ParIndx] = Par->getErrorHi(); + ErrL[ParIndx] = -1*Par->getErrorLo(); + std::cout << Par->GetName() << " " << ParIndx << " " << Res[ParIndx] << " + " << ErrH[ParIndx] << " - " << ErrL[ParIndx]<< std::endl; + } +} + + +void fillVectors(TTree* fr, double Res[][nPars], double ErrH[][nPars], double ErrL[][nPars], + double Diff[][nPars], double DErrH[][nPars], double DErrL[][nPars], + double genRes[nPars], double genErrH[nPars], double genErrL[nPars], + int fileN, int nentry=0){ + + for (int ParIndx=0; ParIndxSetBranchAddress(Form("%s_best",ParName[ParIndx].c_str()),&val); + fr->SetBranchAddress(Form("%s_high",ParName[ParIndx].c_str()),&ci[0]); + fr->SetBranchAddress(Form("%s_low" ,ParName[ParIndx].c_str()),&ci[1]); + + fr->GetEntry(0); + Res[fileN][ParIndx] = val; + ErrH[fileN][ParIndx] = ci[0]-val; + ErrL[fileN][ParIndx] = val - ci[1]; + if (nentry==-1){ + ErrH[fileN][ParIndx] = ci[0]; + ErrL[fileN][ParIndx] = ci[1]; + } + std::cout << ParName[ParIndx].c_str() << " " << ParIndx << " " << val << " + " << ErrH[fileN][ParIndx] << " - " << ErrL[fileN][ParIndx]<< std::endl; + Diff[fileN][ParIndx] = Res[fileN][ParIndx]-genRes[ParIndx]; + DErrH[fileN][ParIndx] = sqrt(ErrH[fileN][ParIndx]*ErrH[fileN][ParIndx]+genErrL[ParIndx]*genErrL[ParIndx]); + DErrL[fileN][ParIndx] = sqrt(ErrL[fileN][ParIndx]*ErrL[fileN][ParIndx]+genErrH[ParIndx]*genErrH[ParIndx]); + } +} + + +void setGraphProperties(TGraphAsymmErrors* Gr, int icolor, string label, string type, int marker){ + + int addc = 0; + if (type.find("WT") != std::string::npos) addc = 1; + if (type.find("CT") != std::string::npos) addc = 2; + + Gr -> SetLineColor(colorlist[icolor] + addc); + Gr -> SetMarkerColor(colorlist[icolor] + addc); + Gr -> SetMarkerStyle(marker + icolor + addc); + Gr -> SetName(Form("Gr%s_%s" , type.c_str(), label.c_str()) ); + Gr -> SetLineWidth(2); +} + +void setAxisProperties(TH1F* auxE2){ + auxE2->GetYaxis()->SetRangeUser(-1*diffMax,diffMax); + auxE2->SetStats(kFALSE); + auxE2->SetLineColor(1); + auxE2->GetXaxis()->SetLabelSize( 0.2); + auxE2->GetXaxis()->SetTickLength(0.1); + auxE2->GetYaxis()->SetTitle("(RECO - GEN)"); + auxE2->GetYaxis()->SetTitleSize(0.07); + auxE2->GetYaxis()->SetTitleOffset(0.42); + auxE2->GetYaxis()->SetLabelSize(0.08); + auxE2->GetXaxis()->SetNdivisions(800); + + for (int ipar=0; ipar < nPars; ipar++){ + auxE2->GetXaxis()->SetBinLabel(ipar+1, Form("%s", ParTitle[ipar].c_str())); + } +} + +void plotFitResultsResonantBin(int parity, bool plotRECO, std::vector files, std::vector labels, bool datalike) +{ + + double pVal [nPars]; + double pErrH [nPars]; + double pErrL [nPars]; + + for (int i=0; i Gr, GrDiff; + + TFile* finGen = TFile::Open(Form("/eos/user/a/aboletti/BdToKstarMuMu/simFitResults/fitResult_genMC_penalty.root")); + if ( !finGen || finGen->IsZombie() ) { + cout<<"Missing gen file: fitResults/fitResult_genMC.root"<Get(Form("ws_b%ip%i_s0_pow1.0",thebin,parity)); + if (wspRes && !wspRes->IsZombie()) { + RooFitResult* fitResultGen = (RooFitResult*)wspRes->obj("fitResult"); + if (fitResultGen && !fitResultGen->IsZombie() && fitResultGen->status()==0 && fitResultGen->covQual()==3) { + fillVectors(fitResultGen, genRes, genErrH, genErrL); + } + } + + // first fill for gen results //roofitresult version +// RooFitResult* fitResultGen = (RooFitResult*)finGen->Get(Form("fitResult_b%ip%i",thebin,parity)); +// if (fitResultGen && !fitResultGen->IsZombie() && fitResultGen->status()==0 && fitResultGen->covQual()==3) { +// fillVectors(fitResultGen, genRes, genErrH, genErrL); +// } + finGen->Close(); + + string label = ""; + for (unsigned int ifile = 0; ifile < n_files; ifile++) { + label.clear(); label.assign(labels[ifile]); + + std::fill(std::begin(Res[ifile]), std::end(Res[ifile]), 0); + + TFile* finFullReco; + std::string filename = files[ifile]; + if (plotRECO) finFullReco = TFile::Open( filename.c_str() ); + + // fill for reco angular fit + if ( plotRECO && finFullReco && !finFullReco->IsZombie() ) { + TTree* ttree ; + ttree = (TTree*)finFullReco->Get("fitResultsTree"); + fillVectors(ttree, Res, ErrH, ErrL, Diff, DiffErrH, DiffErrL, genRes, genErrH, genErrL, ifile, datalike ? 0 : -1); + + if ( plotRECO && finFullReco ) finFullReco->Close(); + } + + // create TGraphs with fit results + Gr .push_back( new TGraphAsymmErrors(nBins, pVal, Res[ifile], pErrH, pErrL, ErrL[ifile], ErrH[ifile]) ); + setGraphProperties(Gr[ifile], ifile, labels[ifile], "Full", 22); + + // create TGraphs with difference to GEN Results + GrDiff.push_back( new TGraphAsymmErrors(nBins, pValReco[ifile], Diff[ifile], pErrHReco, pErrLReco, DiffErrL[ifile], DiffErrH[ifile]) ); + setGraphProperties(GrDiff[ifile], ifile, labels[ifile], "Diff", 22); + + } + + + TGraphAsymmErrors* GrGen = new TGraphAsymmErrors(nPars,pVal,genRes,pErrH,pErrL,genErrL,genErrH); + GrGen->SetName("GrGen"); + + GrGen->SetTitle( Form("Fit result comparison for bin %i ", thebin)); + GrGen->GetYaxis()->SetTitleSize(20); + GrGen->GetYaxis()->SetTitleFont(43); + GrGen->GetYaxis()->SetTitleOffset(1.55); + + GrGen->GetXaxis()->SetRangeUser(0, nPars); + GrGen->GetYaxis()->SetRangeUser(-1, 1); + GrGen->SetLineWidth(2); + + // Legend + TLegend *leg; + leg = new TLegend(0.15,0.1,0.4,0.3); + leg->SetBorderSize(0); + leg->SetFillColor(0); + leg->SetTextSize(0.032); + leg->AddEntry(GrGen,"Fit to generation-level events","lep"); + + // Zero line + TLine *line = new TLine(0,0,nPars,0); + line->SetLineColor(14); + line->SetLineStyle(7); + + c[0] = new TCanvas("c","c",800,800); + c[0]->cd(); + TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0); + pad1->SetBottomMargin(0); + pad1->Draw(); + pad1->cd(); + + GrGen->Draw("AP"); + for (unsigned int ifile = 0; ifile < files.size() ; ifile++) { + Gr[ifile] -> Draw("P"); + leg->AddEntry(Gr[ifile], Form("Fit to reconstructed events, %s", labels[ifile].c_str()), "lep"); + } + + leg->Draw(); + + GrGen->GetYaxis()->SetLabelSize(0.); + TGaxis *axis = new TGaxis( 0 , -1, + 0 , 1, + -1, 1, + 510, ""); + axis->SetName("leg1"); + axis->SetLabelFont(43); + axis->SetLabelSize(15); + axis->Draw(); + + + // plot difference wrt GEN results + c[0]->cd(); + TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3); + pad2->SetTopMargin(0); + pad2->SetBottomMargin(0.25); + pad2->Draw(); + pad2->cd(); + + // first create axis + TH1F* auxE2 = new TH1F("auxE2", "", nPars, 0, nPars); + setAxisProperties(auxE2); + auxE2->Draw(); + + for (unsigned int ifile = 0; ifile < files.size() ; ifile++) { + if (plotRECO) GrDiff[ifile] -> Draw("P"); + } + line->Draw(); + + string confString = "plotSimFit_d/fitResult_"; + if (plotRECO) confString = confString + "recoRes"; + c[0]->SaveAs( (confString + stat + Form("_bin%d_4D_XGBv8.pdf",thebin)).c_str() ); +} + +void compareSimFitResultsResonantBins(int parity, bool plotRECO = true, bool datalike = false) +{ + + gROOT->SetBatch(true); + if ( parity<0 || parity>1 ) return; + + std::vector files, labels; + //MC stat 3D +// files.push_back(Form("simFitResults/xgbv8/simFitResult_recoMC_fullAngular201620172018_MCStat_b%dp1_XGBv8.root",thebin)); +// labels.push_back("3D-fit"); +// files.push_back(Form("simFitResults/xgbv8/simFitResult_recoMC_fullAngular2016_MCStat_b%dp1_XGBv8.root",thebin)); +// labels.push_back("3D-fit, 2016"); +// files.push_back(Form("simFitResults/xgbv8/simFitResult_recoMC_fullAngular2017_MCStat_b%dp1_XGBv8.root",thebin)); +// labels.push_back("3D-fit, 2017"); +// files.push_back(Form("simFitResults/xgbv8/simFitResult_recoMC_fullAngular2018_MCStat_b%dp1_XGBv8.root",thebin)); +// labels.push_back("3D-fit, 2018"); + + //MC stat 4D + gStyle->SetPalette(77); + files.push_back(Form("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass201620172018_MCStat_b%dp1_XGBv8.root",thebin)); + labels.push_back("4D-fit, simultaneous 3 years"); + files.push_back(Form("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass2016_MCStat_b%d_XGBv8.root",thebin)); + labels.push_back("4D-fit, 2016"); + files.push_back(Form("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass2017_MCStat_b%d_XGBv8.root",thebin)); + labels.push_back("4D-fit, 2017"); + files.push_back(Form("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass2018_MCStat_b%d_XGBv8.root",thebin)); + labels.push_back("4D-fit, 2018"); + + // data stat, subsample 2 +// gStyle->SetPalette(77); +// files.push_back("simFitResults/newphi/simFitResult_recoMC_fullAngular201620172018_dataStat_newPhi_b1.root"); +// labels.push_back("3D-fit"); +// files.push_back("simFitResults4d/simFitResult_recoMC_fullAngularMass201620172018_dataStat_b0.root"); +// labels.push_back("4D-fit"); +// files.push_back("simFitResults4d/simFitResult_recoMC_fullAngularMass_toybkg201620172018_dataStat-1_update_b0.root"); +// labels.push_back("toybkg"); + +// // data stat, 100 subsamples +// gStyle->SetPalette(71); +// files.push_back("simFitResults/simFitResult_recoMC_fullAngular201620172018_dataStat-*_b0.root"); +// labels.push_back("3D-fit"); +// files.push_back("simFitResults/simFitResult_recoMC_fullAngular201620172018_dataStat-*_b0.root"); +// labels.push_back("4D-fit"); +// files.push_back("simFitResults4d/simFitResult_recoMC_fullAngularMass_toybkg201620172018_dataStat-*_b0.root"); +// labels.push_back("toybkg"); + + if (datalike) diffMax = 0.3999; + + colorlist[0] = TColor::GetColorPalette(80); + colorlist[1] = TColor::GetColorPalette(160); + colorlist[2] = TColor::GetColorPalette(240); + colorlist[3] = TColor::GetColorPalette(10); + colorlist[4] = TColor::GetColorPalette(200); + + plotFitResultsResonantBin(parity, plotRECO, files, labels, datalike); + +} diff --git a/computeCoverage_dataFitToys.cc b/computeCoverage_dataFitToys.cc new file mode 100644 index 0000000..971155c --- /dev/null +++ b/computeCoverage_dataFitToys.cc @@ -0,0 +1,185 @@ +#include +#include +#include +#include +#include +#include + +using namespace std; + +string label = "test"; + +static const int nPars = 8; +string parName [nPars] = {"Fl","P1","P2","P3","P4p","P5p","P6p","P8p"}; +string parTitle[nPars] = {"F_{L}","P_{1}","P_{2}","P_{3}","P'_{4}","P'_{5}","P'_{6}","P'_{8}"}; + +void computeCoverage_dataFitToys (int q2Bin = 2) +{ + + int parity = 1; + + vector< vector > vCover (nPars); + vector< vector > vCoErr (nPars); + + vector< vector > vX (nPars); + vector< vector > vXe (nPars); + + string filename_data = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root", q2Bin); + auto file_data = new TFile(filename_data.c_str(),"READ"); + if (!file_data) { + cout<<"File not found: "<Get("fitResultsTree"); + + double chi2 = 0; + double chi = 0; + int ndf = 0; + + TChain MINOS_output("MINOS_output",""); + string filename = Form("simFitResults4d/simFitResult_toy_fullAngularMass_Swave_201620172018_s*_b%i-XGBv8_unbl4.root/fitResultsTree",q2Bin); + // filename = "/afs/cern.ch/work/d/dini/public/Coverage/" + filename; + MINOS_output.Add(filename.c_str()); + + int ntoys = MINOS_output.GetEntries(); + cout<<"Toys: "< vRef (nPars); + vector vRes (nPars); + vector vHigh(nPars); + vector vLow (nPars); + vector nIn (nPars); + vector hdist(0); + vector hdout(0); + for (int iPar=0; iParSetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vRef[iPar]); + + MINOS_output.SetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vRes[iPar]); + MINOS_output.SetBranchAddress(Form("%s_high",parName[iPar].c_str()),&vHigh[iPar]); + MINOS_output.SetBranchAddress(Form("%s_low" ,parName[iPar].c_str()),&vLow [iPar]); + + hdist.push_back( new TH1F(Form("hdist%i",iPar),Form(";%s",parTitle[iPar].c_str()),100,-1,1) ); + hdout.push_back( new TH1F(Form("hdout%i",iPar),Form(";%s",parTitle[iPar].c_str()),100,-1,1) ); + nIn[iPar] = 0; + + } + + result_data->GetEntry(0); + + int testPar = 0; + TLine line2 (vRef[testPar],0,vRef[testPar],300); + // vRef[7]=0.125063; + // vRef[2]=0.0895366; + // vRef[5]=-0.360472; + // vRef[6]=0.0887163; + // vRef[7]=0.137497; + // vRef[0]=0.706922; + vector vTestX (0); + vector vTestXeh (0); + vector vTestXel (0); + vector vTestY (0); + vector vTestYeh (0); + vector vTestYel (0); + vector voTestX (0); + vector voTestXeh (0); + vector voTestXel (0); + vector voTestY (0); + vector voTestYeh (0); + vector voTestYel (0); + + for (int iEn=0; iEn vRes[testPar]+vLow[testPar] ) { + // vTestX.push_back(iEn+1); + // vTestXeh.push_back(0); + // vTestXel.push_back(0); + // vTestY.push_back(vRes[testPar]); + // vTestYeh.push_back(vHigh[testPar]); + // vTestYel.push_back(vLow[testPar]); + // } else { + // voTestX.push_back(iEn+1); + // voTestXeh.push_back(0); + // voTestXel.push_back(0); + // voTestY.push_back(vRes[testPar]); + // voTestYeh.push_back(vHigh[testPar]); + // voTestYel.push_back(vLow[testPar]); + // } + for (int iPar=0; iParFill(vRes[iPar]); + if ( vRef[iPar] < vRes[iPar]+vHigh[iPar] && vRef[iPar] > vRes[iPar]+vLow[iPar] ) { + hdout[iPar]->Fill(vRes[iPar]); + ++nIn[iPar]; + } + } + } + + for (int iPar=0; iPar grCover (nPars); + for (int iPar=0; iParSetName(Form("grCov%s",parName[iPar].c_str())); + grCover[iPar]->SetMarkerStyle(20); + } + + gStyle->SetOptStat(0); + cCover.cd(1)->SetTicky(2); + string plotTitle = Form("Toy coverage study - q2 bin %i;;coverage",q2Bin); + auto hLab = new TH1S ("hLab",plotTitle.c_str(),nPars,0,nPars); + for (int iBin=0; iBinGetXaxis()->SetBinLabel(iBin+1,parTitle[iBin].c_str()); + hLab->SetMinimum(38); + hLab->SetMaximum(88); + hLab->GetYaxis()->SetTickLength(0.006); + hLab->Draw(); + + TLine line (0,68.27,nPars,68.27); + line.SetLineStyle(2); + line.SetLineColor(15); + line.Draw(); + + for (int iPar=0; iParDraw("P"); + cout<GetMean()<<" "<GetRMS()<GetXaxis()->SetRangeUser(hdist[testPar]->GetMean()-5*hdist[testPar]->GetRMS(),hdist[testPar]->GetMean()+5*hdist[testPar]->GetRMS()); + // hdist[testPar]->Draw(); + // hdout[testPar]->SetLineColor(8); + // hdout[testPar]->Draw("same"); + // line2.SetLineColor(2); + // line2.Draw(); + + cCover.SaveAs(Form("plotSimFit4d_d/coverage_datatoy_b%i.pdf",q2Bin)); + +} diff --git a/computeCoverage_dataFitToys_allbins.cc b/computeCoverage_dataFitToys_allbins.cc new file mode 100644 index 0000000..bca9a60 --- /dev/null +++ b/computeCoverage_dataFitToys_allbins.cc @@ -0,0 +1,192 @@ +#include +#include +#include +#include +#include +#include + +using namespace std; + +string label = "test"; + +static const int nPars = 8; +string parName [nPars] = {"Fl","P1","P2","P3","P4p","P5p","P6p","P8p"}; +string parTitle[nPars] = {"F_{L}","P_{1}","P_{2}","P_{3}","P'_{4}","P'_{5}","P'_{6}","P'_{8}"}; + +vector binNumb = {0,1,2,3,5,7}; + +void computeCoverage_dataFitToys_allbins () +{ + + vector< vector > vCover (nPars); + vector< vector > vCoErr (nPars); + + vector< vector > vX (nPars); + vector< vector > vXe (nPars); + + double chi2 = 0; + double chi = 0; + int ndf = 0; + + for (int iBin = 0; iBinGet("fitResultsTree"); + + TChain MINOS_output("MINOS_output",""); + string filename = Form("simFitResults4d/simFitResult_toy_fullAngularMass_Swave_201620172018_s*_b%i-XGBv8_unbl4.root/fitResultsTree",q2Bin); + // filename = "/afs/cern.ch/work/d/dini/public/Coverage/" + filename; + MINOS_output.Add(filename.c_str()); + + int ntoys = MINOS_output.GetEntries(); + cout<<"Toys: "< vRef (nPars); + vector vRes (nPars); + vector vHigh(nPars); + vector vLow (nPars); + vector nIn (nPars); + // vector hdist(0); + // vector hdout(0); + for (int iPar=0; iParSetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vRef[iPar]); + + MINOS_output.SetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vRes[iPar]); + MINOS_output.SetBranchAddress(Form("%s_high",parName[iPar].c_str()),&vHigh[iPar]); + MINOS_output.SetBranchAddress(Form("%s_low" ,parName[iPar].c_str()),&vLow [iPar]); + + // hdist.push_back( new TH1F(Form("hdist%i",iPar),Form(";%s",parTitle[iPar].c_str()),100,-1,1) ); + // hdout.push_back( new TH1F(Form("hdout%i",iPar),Form(";%s",parTitle[iPar].c_str()),100,-1,1) ); + nIn[iPar] = 0; + + } + + result_data->GetEntry(0); + + int testPar = 0; + TLine line2 (vRef[testPar],0,vRef[testPar],300); + // vRef[7]=0.125063; + // vRef[2]=0.0895366; + // vRef[5]=-0.360472; + // vRef[6]=0.0887163; + // vRef[7]=0.137497; + // vRef[0]=0.706922; + vector vTestX (0); + vector vTestXeh (0); + vector vTestXel (0); + vector vTestY (0); + vector vTestYeh (0); + vector vTestYel (0); + vector voTestX (0); + vector voTestXeh (0); + vector voTestXel (0); + vector voTestY (0); + vector voTestYeh (0); + vector voTestYel (0); + + for (int iEn=0; iEn vRes[testPar]+vLow[testPar] ) { + // vTestX.push_back(iEn+1); + // vTestXeh.push_back(0); + // vTestXel.push_back(0); + // vTestY.push_back(vRes[testPar]); + // vTestYeh.push_back(vHigh[testPar]); + // vTestYel.push_back(vLow[testPar]); + // } else { + // voTestX.push_back(iEn+1); + // voTestXeh.push_back(0); + // voTestXel.push_back(0); + // voTestY.push_back(vRes[testPar]); + // voTestYeh.push_back(vHigh[testPar]); + // voTestYel.push_back(vLow[testPar]); + // } + for (int iPar=0; iParFill(vRes[iPar]); + if ( vRef[iPar] < vRes[iPar]+vHigh[iPar] && vRef[iPar] > vRes[iPar]+vLow[iPar] ) { + // hdout[iPar]->Fill(vRes[iPar]); + ++nIn[iPar]; + } + } + } + + for (int iPar=0; iPar grCover (nPars); + for (int iPar=0; iParSetName(Form("grCov%s",parName[iPar].c_str())); + grCover[iPar]->SetMarkerStyle(20); + } + + int nBins = binNumb.size(); + gStyle->SetOptStat(0); + cCover.cd(1)->SetTicky(2); + string plotTitle = "Toy coverage study;;coverage"; + auto hLab = new TH1S ("hLab",plotTitle.c_str(),nBins,0,nBins); + for (int iBin=0; iBinGetXaxis()->SetBinLabel(iBin+1,Form("q2-bin %i",binNumb[iBin])); + hLab->SetMinimum(48); + hLab->SetMaximum(88); + hLab->GetYaxis()->SetTickLength(0.006); + hLab->Draw(); + + TLine line (0,68.27,nBins,68.27); + line.SetLineStyle(2); + line.SetLineColor(15); + line.Draw(); + + for (int iPar=0; iParDraw("P"); + // cout<GetMean()<<" "<GetRMS()<GetXaxis()->SetRangeUser(hdist[testPar]->GetMean()-5*hdist[testPar]->GetRMS(),hdist[testPar]->GetMean()+5*hdist[testPar]->GetRMS()); + // hdist[testPar]->Draw(); + // hdout[testPar]->SetLineColor(8); + // hdout[testPar]->Draw("same"); + // line2.SetLineColor(2); + // line2.Draw(); + + cCover.SaveAs("plotSimFit4d_d/coverage_datatoy.pdf"); + +} diff --git a/computeSystematicsFromDataFit.cc b/computeSystematicsFromDataFit.cc new file mode 100644 index 0000000..d728bbb --- /dev/null +++ b/computeSystematicsFromDataFit.cc @@ -0,0 +1,206 @@ +#include +#include +#include +#include +#include +#include + +using namespace std; + +static const int nBins = 8; +float binBorders [nBins+1] = { 1.1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; + +static const int nUncBins = 100; +double binsUnc [nUncBins+1]; + +static const int nPars = 8; +string parName [nPars] = {"Fl","P1","P2","P3","P4p","P5p","P6p","P8p"}; +string parTitle[nPars] = {"F_{L}","P_{1}","P_{2}","P_{3}","P'_{4}","P'_{5}","P'_{6}","P'_{8}"}; +double parMin [nPars] = {0,-1,-0.5,-0.5,-1*sqrt(2),-1*sqrt(2),-1*sqrt(2),-1*sqrt(2)}; +double parMax [nPars] = {1, 1, 0.5, 0.5, sqrt(2), sqrt(2), sqrt(2), sqrt(2)}; + +static const int nQuant = 4; +double quantPerc [nQuant] = {0.025,0.16,0.84,0.975}; + +int colors [12] = { 633, 417, 879, 857, 839, 801, 921, 607, 807, 419, 907, 402 }; +// int colors [13] = { 633, 417, 879, 857, 839, 887, 801, 921, 607, 807, 419, 907, 402 }; + +static const int nVariations = 1; +int aq2stat [nVariations] = { 0 }; + +double diffMax = 0.0499; + +void computeSystematicsFromDataFit () +{ + + gROOT->SetBatch(true); + gStyle->SetOptStat(0); + + vector vq2Bins (0); + vector table (0); + + binsUnc[0]=0.006; + for (int i=0; i vNominal(nPars); + vector vHigh(nPars); + vector vLow (nPars); + + vector vSyst(nPars); + vector vSystH(nPars); + vector vSystL(nPars); + + if (q2Bin==4 || q2Bin==6) { + continue; + } + + vector vBias(nPars); + + int q2stat = aq2stat[0]; + vq2Bins.push_back(q2stat); + + + TChain fitResultsTree ("fitResultsTree",""); + // mass shape syst + // string filename = Form("/eos/user/a/aboletti/BdToKstarMuMu/massShapeSyst/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_f1.06_unbl4.root",q2Bin); + // mistag syst + string filename = Form("simFitResults4d/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4_mistagSyst_exp.root",q2Bin); + if (q2Bin==7){ + filename = "/afs/cern.ch/work/d/dini/public/LastSystStudies/Bin7mfracFree/simFitResults4d/simFitResult_data_fullAngularMass_Swave_201620172018_b7-XGBv8_unbl4_mistagSyst_exp.root"; + } + fitResultsTree.Add(filename.c_str()); + + string filename_nominal = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); +// string filename_nominal = Form("simFitResults4d/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_wp90_unbl4.root",q2Bin); + TFile* filein_nominal = TFile::Open(filename_nominal.c_str()); + TTree* fitResultsTree_nominal = (TTree*)filein_nominal->Get("fitResultsTree"); + if (!fitResultsTree_nominal || fitResultsTree_nominal->GetEntries() != 1) { + cout<<"Error, unexpected number of entries in fitResultsTree in file: "<SetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vNominal[iPar]); + fitResultsTree_nominal->SetBranchAddress(Form("%s_high",parName[iPar].c_str()),&vHigh[iPar]); + fitResultsTree_nominal->SetBranchAddress(Form("%s_low" ,parName[iPar].c_str()),&vLow [iPar]); + + fitResultsTree.SetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vSyst[iPar]); + fitResultsTree.SetBranchAddress(Form("%s_high",parName[iPar].c_str()),&vSystH[iPar]); + fitResultsTree.SetBranchAddress(Form("%s_low" ,parName[iPar].c_str()),&vSystL[iPar]); + } + + fitResultsTree_nominal->GetEntry(0); + fitResultsTree.GetEntry(0); +// for (int iPar=0; iPar0.0005) + line = line + " & 0.001"; + else + line = line + Form(" & %.3f",(vBias[iPar])); + } + if (table.size()==0) { + firstLine = firstLine + " \\\\ \\hline"; + table.push_back(firstLine); + } + line = line + " \\\\"; + table.push_back(line); + + // Plot P-wave parameters + + double x[nPars]; + double xe[nPars]; + double xNom[nPars]; + double xeNom[nPars]; + + double yNominal [nPars]; + double ylNominal[nPars]; + double yhNominal[nPars]; + double y [nPars]; + double yl[nPars]; + double yh[nPars]; + + for (int iPar=0; iParGetXaxis()->SetBinLabel(iBin+1,parName[iBin].c_str()); + double yRange = 0.25; + hLab->SetMinimum(-1*yRange); + hLab->SetMaximum(yRange); + hLab->Draw(); + + auto grNominal = new TGraphAsymmErrors(nPars,xNom,yNominal,xeNom,xeNom,ylNominal,yhNominal); + grNominal->SetName("grNominal"); + grNominal->SetFillColor(429); + grNominal->Draw("2"); + + TLine linex (0,0,nPars,0); + linex.SetLineWidth(2); + linex.SetLineStyle(2); + linex.SetLineColor(13); + linex.Draw(); + + TGraphAsymmErrors* gr= new TGraphAsymmErrors(nPars,xNom,y,xeNom,xeNom,yl,yh); + gr->SetLineColor(1); + gr->SetLineWidth(2); + gr->Draw("P"); + + TLegend leg (0.15,0.75,0.4,0.9); + leg.AddEntry(grNominal,"Nominal fit","f"); + leg.AddEntry(gr,"sideband stat syst","lep"); +// leg.AddEntry(gr,"mass shape syst","lep"); +// leg.AddEntry(gr,"mistag fraction syst","lep"); + leg.Draw(); + +// canv.SaveAs(Form("plotSyst/simfit_syst_mistag_results_xgbv8_wp90_bin%d.pdf", q2Bin)); +// canv.SaveAs(Form("plotSyst/simfit_syst_sidebandStat_results_xgbv8_wp90_bin%d.pdf", q2Bin)); + + } // end n bins + + cout<<"===== Formatted bias table ========"<> swapCut = { + {6, {0.13, 0.15, 0.07, 0.43, -0.45, 2.0}}, + {7, {0.07, 0.13, 0.11, 0.45, -0.65, 1.9}}, + {8, {0.09, 0.13, 0.07, 0.23, -0.45, 1.0}} }; TCanvas* c [nBins]; -void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool plot = false) +void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 8, bool plot = false) { // year format: [6] for 2016 // [7] for 2017 @@ -24,7 +31,7 @@ void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool pl if ( year<6 || year>8 ) return; string XGBstr = ""; - if (XGBv>0 && XGBv<6) XGBstr = Form("_XGBv%i",XGBv); + if (!data && (XGBv==8 || (XGBv>0 && XGBv<6))) XGBstr = Form("_XGBv%i",XGBv); bool isJpsi = false; bool isPsi = false; @@ -56,6 +63,7 @@ void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool pl for (int i=0; i0?Form("-XGBv%i",XGBv):""):"data",year)); auto t_num = (TTree*)f_num->Get("ntuple"); - if (data==0 && XGBv>9) { - t_num->AddFriend("wTree",Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/MC-%s-TMVAv%i/201%i.root",isJpsi?"Jpsi":(isPsi?"Psi":"LMNR"),XGBv-10,year)); - XGBstr = Form("_TMVAv%i",XGBv-10); - } else if (data==0 && XGBv>5) { - t_num->AddFriend("wTree",Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/MC-%s-XGBv%i/201%i.root",isJpsi?"Jpsi":(isPsi?"Psi":"LMNR"),XGBv,year)); - XGBstr = Form("_XGBv%i",XGBv); - } + // if (data==0 && XGBv>9) { + // t_num->AddFriend("wTree",Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/MC-%s-TMVAv%i/201%i.root",isJpsi?"Jpsi":(isPsi?"Psi":"LMNR"),XGBv-10,year)); + // XGBstr = Form("_TMVAv%i",XGBv-10); + // } else if (data==0 && XGBv>5) { + // t_num->AddFriend("wTree",Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/MC-%s-XGBv%i/201%i.root",isJpsi?"Jpsi":(isPsi?"Psi":"LMNR"),XGBv,year)); + // XGBstr = Form("_XGBv%i",XGBv); + // } // TChain* t_num = new TChain(); // if (data==0 && isLMNR) @@ -118,6 +126,58 @@ void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool pl t_num->SetBranchAddress( "passB0Psi_jpsi", &passB0Psi_jpsi ); t_num->SetBranchAddress( "passB0Psi_psip", &passB0Psi_psip ); + // Anti-swap cut + double kstTrkmPt = 0; + double kstTrkmEta = 0; + double kstTrkmPhi = 0; + double kstTrkpPt = 0; + double kstTrkpEta = 0; + double kstTrkpPhi = 0; + double mumPt = 0; + double mumEta = 0; + double mumPhi = 0; + double mupPt = 0; + double mupEta = 0; + double mupPhi = 0; + double selected_swapped_b0_mass = 0; + double selected_swapped_kstar_mass = 0; + double selected_swapped_mumu_mass = 0; + Long64_t charge_pf_swapped_track = 0; + + t_num->SetBranchStatus("kstTrkmPt",1); + t_num->SetBranchStatus("kstTrkmEta",1); + t_num->SetBranchStatus("kstTrkmPhi",1); + t_num->SetBranchStatus("kstTrkpPt",1); + t_num->SetBranchStatus("kstTrkpEta",1); + t_num->SetBranchStatus("kstTrkpPhi",1); + t_num->SetBranchStatus("mumPt",1); + t_num->SetBranchStatus("mumEta",1); + t_num->SetBranchStatus("mumPhi",1); + t_num->SetBranchStatus("mupPt",1); + t_num->SetBranchStatus("mupEta",1); + t_num->SetBranchStatus("mupPhi",1); + t_num->SetBranchStatus("selected_swapped_b0_mass",1); + t_num->SetBranchStatus("selected_swapped_kstar_mass",1); + t_num->SetBranchStatus("selected_swapped_mumu_mass",1); + t_num->SetBranchStatus("charge_pf_swapped_track",1); + + t_num->SetBranchAddress("kstTrkmPt" , &kstTrkmPt ); + t_num->SetBranchAddress("kstTrkmEta" , &kstTrkmEta ); + t_num->SetBranchAddress("kstTrkmPhi" , &kstTrkmPhi ); + t_num->SetBranchAddress("kstTrkpPt" , &kstTrkpPt ); + t_num->SetBranchAddress("kstTrkpEta" , &kstTrkpEta ); + t_num->SetBranchAddress("kstTrkpPhi" , &kstTrkpPhi ); + t_num->SetBranchAddress("mumPt" , &mumPt ); + t_num->SetBranchAddress("mumEta" , &mumEta ); + t_num->SetBranchAddress("mumPhi" , &mumPhi ); + t_num->SetBranchAddress("mupPt" , &mupPt ); + t_num->SetBranchAddress("mupEta" , &mupEta ); + t_num->SetBranchAddress("mupPhi" , &mupPhi ); + t_num->SetBranchAddress("selected_swapped_b0_mass" , &selected_swapped_b0_mass ); + t_num->SetBranchAddress("selected_swapped_kstar_mass", &selected_swapped_kstar_mass); + t_num->SetBranchAddress("selected_swapped_mumu_mass" , &selected_swapped_mumu_mass ); + t_num->SetBranchAddress("charge_pf_swapped_track" , &charge_pf_swapped_track ); + // B0-kinematic variables // double recoB0pT, recoB0eta; // t_num->SetBranchAddress( "bPt" , &recoB0pT ); @@ -155,7 +215,11 @@ void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool pl // MC weights double XGBweight = 1; float fXGBweight = 1; - if (XGBv>9) { + if (XGBv==8) { + t_num->SetBranchStatus("XGBv8",1); + t_num->SetBranchAddress( "XGBv8", &fXGBweight ); + } + else if (XGBv>9) { t_num->SetBranchStatus("MCw",1); t_num->SetBranchAddress( "MCw", &XGBweight ); } @@ -208,7 +272,18 @@ void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool pl if (xBin<0) continue; // apply cut for bin 4 - if (isJpsi && XCut>0) continue; + if (isJpsi && XCut>0) continue; + + // anti-swap cut + if (fabs(selected_swapped_mumu_mass-PDGJpsiMass)swapCut[year][4] && + (mumPt-kstTrkmPt)/kstTrkmPt0 && (deltaR(kstTrkpEta,kstTrkpPhi,mupEta,mupPhi)swapCut[year][4] && + (mupPt-kstTrkpPt)/kstTrkpPt 1.0*counter*numEntries/100 ) { @@ -414,6 +489,17 @@ void createDataset(int year, int q2Bin = -1, int data = 0, int XGBv = 4, bool pl // apply cut for bin 4 if (isJpsi && XCut>0) continue; + // anti-swap cut + if (fabs(selected_swapped_mumu_mass-PDGJpsiMass)swapCut[year][4] && + (mumPt-kstTrkmPt)/kstTrkmPt0 && (deltaR(kstTrkpEta,kstTrkpPhi,mupEta,mupPhi)swapCut[year][4] && + (mupPt-kstTrkpPt)/kstTrkpPt 1.0*counter*numEntries/100 ) { cout<TMath::Pi()) { + if (deltaPhi>0) deltaPhi = deltaPhi - 2*TMath::Pi(); + else deltaPhi = deltaPhi + 2*TMath::Pi(); + } + float deltaEta = eta1 - eta2; + return sqrt( deltaEta*deltaEta + deltaPhi*deltaPhi ); +} diff --git a/generateToy.cc b/generateToy.cc new file mode 100644 index 0000000..4fa84a4 --- /dev/null +++ b/generateToy.cc @@ -0,0 +1,104 @@ +#include +#include + +#include +#include +#include +#include +#include + +#include "utils.h" +#include "PdfSigRTMass.h" +#include "PdfSigWTMass.h" +#include "PdfSigAngMass.h" +#include "ShapeSigAng.h" + +#include "BoundCheck.h" +#include "BoundDist.h" +#include "Penalty.h" +#include "Fitter.h" +#include "RooBernsteinSideband.h" + +using namespace RooFit; +using namespace std; + +vector parNames = {"Fl", "P1", "P2", "P3", "P4p", "P5p", "P6p", "P8p"}; + +int main(int argc, char** argv) +{ + + int q2Bin = 2; + int parity = 0; + + if ( argc > 1 ) q2Bin = atoi(argv[1]); + if ( argc > 2 ) parity = atoi(argv[2]); + + string finname = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); + auto fin = new TFile(finname.c_str(),"READ"); + if (!fin) { + cout<<"File not found: "<Get("wsp_out"); + if (!wsp) { + cout<<"Workspace not found in file: "<Get("fitResultsTree"); + if (!tin) { + cout<<"TTree not found in file: "< dataResults(parNames.size(),0); + for (uint iPar=0; iParSetBranchAddress((parNames[iPar]+"_best").c_str(),&dataResults[iPar]); + tin->GetEntry(0); + + cout<<"From tree"; + for (uint iPar=0; iPardata("allcombData"); + auto protopdf = (RooSimultaneous*)wsp->pdf("simPdf"); + auto observables = protopdf->getObservables(protodata); + auto protoparameters = protopdf->getParameters(observables); + + cout<<"From PDF "; + for (uint iPar=0; iParfind(parNames[iPar].c_str()))->getVal(); + cout<reduce(*observables,Form("sample==sample::data%i_subs0",year))->numEntries(); + + auto singlepdf = protopdf->getPdf(Form("data%i_subs0",year)); + + cout<<"Generating "<generate(*observables, nEv, Name(Form("data_b%i",q2Bin)) ); + subTime.Stop(); + cout<<"CPU time: "<Print(); + + auto fout = new TFile( Form( "recoTOYDataset_b%i_%i.root", q2Bin, year), "RECREATE" ); + auto wsp2 = new RooWorkspace(Form("ws_b%ip0",q2Bin),""); + wsp2->import( *toydata ); + wsp2->Write(); + fout->Close(); + + delete wsp2; + delete toydata; + + } + + // singlepdf->fitTo(*toydata); + + // cout<<"From fit "; + // for (uint iPar=0; iPargetParameters(observables))->find(parNames[iPar].c_str()))->getVal(); + // cout< vConfInterLow; /* Custom MINOS error */ std::vector vConfInterHigh; /* Custom MINOS error */ + std::vector vSResult; + std::vector vSErrLow; + std::vector vSErrHigh; + std::vector vFreeSResult; + std::vector vFreeSErrLow; + std::vector vFreeSErrHigh; + Fitter() {} ; Fitter(const char *_name, const char *_title, RooArgList _angPars, @@ -125,6 +137,27 @@ class Fitter { BoundDist* _bound_dist, RooArgSet* _constrVars = nullptr ); + Fitter(const char *_name, const char *_title, + RooArgList _angPars, + RooArgList _sPars, + RooDataSet* _combData, + RooAbsPdf* _simPdf, + RooAbsPdf* _simPdf_penalty, + BoundCheck* _boundary, + BoundDist* _bound_dist, + Penalty* _penTerm, + RooArgSet* _constrVars = nullptr + ); + Fitter(const char *_name, const char *_title, + RooArgList _angPars, + RooArgList _sPars, + RooDataSet* _combData, + RooAbsPdf* _simPdf, + RooAbsPdf* _simPdf_penalty, + BoundCheck* _boundary, + BoundDist* _bound_dist, + RooArgSet* _constrVars = nullptr + ); Fitter(const Fitter& other, const char* name=0) ; virtual Fitter* clone(const char* newname) const { return new Fitter(*this,newname); } inline virtual ~Fitter() { } @@ -145,6 +178,10 @@ class Fitter { void setUnbl(Int_t _unbl = 0) { unbl=_unbl; }; + void PrintMinosScan(Bool_t flag = true) { printMinosScan=flag; }; + + void DisableConstOpt(Bool_t flag = true) { disableConstOpt=flag; }; + ClassDef(Fitter,1) // Code to run the fit and statistical uncertainty }; diff --git a/interface/utils.h b/interface/utils.h index d38f915..93f3592 100644 --- a/interface/utils.h +++ b/interface/utils.h @@ -1,9 +1,11 @@ #include +#include #include #include #include #include #include +#include #include using namespace std; @@ -48,6 +50,17 @@ std::map> constr_width_fac = { {2018, {1.87, 4.17, 1.16, 2.72, 1.27, 2.99, 1.08, 2.13, 1.00, 1.00, 1.78, 1.70, 0.00, 0.00, 1.84, 2.62}} }; + +std::map> deltamv = { + {2016, {0.20,0.10,0.08,0.11}}, + {2017, {0.19,0.09,0.07,0.08}}, + {2018, {0.19,0.09,0.07,0.08}}, +}; + +float binBorders [10] = { 1.1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16, 19}; + + + RooPlot* prepareFrame(RooPlot* frame){ frame->GetYaxis()->SetTitleOffset(1.8); frame->SetMaximum(frame->GetMaximum()*1.15); @@ -198,7 +211,7 @@ std::vector createDatasetInData(int nSample, uint firstSample, uint ("data_"+shortString + Form("_subs%i", is)).c_str(), RooArgSet(vars)); - isample->append(*((RooDataSet*)data->reduce(EventRange(is*subSampStat+1, + isample->append(*((RooDataSet*)data->reduce(EventRange(is*subSampStat, (is+1)*subSampStat)))); datasample.push_back (isample); } @@ -216,3 +229,72 @@ std::vector createDatasetInData(int nSample, uint firstSample, uint return datasample; } +std::vector createDatasetInToy(int nSample, int totSamples, RooWorkspace *ws, + int q2Bin, RooArgSet vars, std::string shortString ) +{ + + if (nSample>=totSamples) + return std::vector(0); + + std::vector datasample; + + RooDataSet* data = (RooDataSet*)ws->data(Form("data_b%i",q2Bin)) ; + + int subSampStat = data->numEntries()/totSamples; + + RooDataSet* isample = new RooDataSet(("data_"+shortString + Form("_subs%i", nSample+1)).c_str(), + ("data_"+shortString + Form("_subs%i", nSample+1)).c_str(), + RooArgSet(vars)); + + isample->append(*((RooDataSet*)data->reduce(EventRange(nSample*subSampStat, + (nSample+1)*subSampStat)))); + datasample.push_back (isample); + + return datasample; +} + +RooGenericPdf* createPdfCuts(int q2Bin, int year, RooRealVar *var, RooRealVar *slo){ + if((q2Bin!=3&&q2Bin!=5&&q2Bin!=7) || (year<2016||year>2018)){ + std::cout< PdfCut not defined for q2Bin=%i, year=%i",q2Bin,year)<%.5f)*(@0-%.5f)/%.5f",m1[i],m2[i],m2[i],binhigh[i]-binlow[i]); + formula[i] = formula[i] + Form(" + (@0<%.5f)*(@0>%.5f)*(%.5f-@0)/%.5f",m3[i],m4[i],m3[i],binhigh[i]-binlow[i]); + } +// + if(q2Bin==3) expression = Form("(@0<%.5f)+(@0>%.5f)",m4[0],m1[0])+formula[0]; + if(q2Bin==5) expression = Form("(@0<%.5f)+(@0>%.5f)*(@0<%.5f)+(@0>%.5f)",m4[1],m1[1],m4[2],m1[2])+formula[1]+formula[2]; + if(q2Bin==7) expression = Form("(@0<%.5f)+(@0>%.5f)",m4[3],m1[3])+formula[3]; + expression="("+expression+")*exp(@0*@1)"; + std::cout< %f, %f, %f, %f",year,deltam[0],deltam[1],deltam[2],deltam[3])<Print(); + return PdfCut; +} diff --git a/plotMassPulls.cc b/plotMassPulls.cc new file mode 100644 index 0000000..377c0eb --- /dev/null +++ b/plotMassPulls.cc @@ -0,0 +1,438 @@ +#include "interface/utils.h" + +string year[3] = {"2016","2017","2018"}; + +void plotMassPulls(int q2Bin = 0, int q2Stat = -1) +{ + gROOT->SetBatch(true); + gStyle->SetOptStat(0); + + bool verbose = true; + + // for data unblinded should be p0, otherwise p1 + string shortString = Form("b%ip0",q2Bin); + string statString = q2Stat<0 ? "" : Form("_b%istat-0",q2Stat); + + static const int nPars = 8; + string parName[nPars] = {"Fl","P1","P2","P3","P4p","P5p","P6p","P8p"}; + + string massParName[] = {"mean_{RT}^{%s}", + "fsig_"+shortString+"_%s"}; + // "slope^{%s}"}; + static const int nMassPars = sizeof(massParName)/sizeof(string); + + std::cout << "bin: " << q2Bin << std::endl; + std::vector massConstParName; + + if (q2Bin >= 5 && q2Bin < 7) + massConstParName = { "f_{M}^{%s}", + "#alpha_{RT1}^{%s}", + "#alpha_{RT2}^{%s}", + "#sigma_{RT1}^{%s}", + "#sigma_{RT2}^{%s}", + "f^{RT%s}", + "n_{RT1}^{%s}", + "n_{RT2}^{%s}", + "#alpha_{WT1}^{%s}", + "#alpha_{WT2}^{%s}", + "#sigma_{WT1}^{%s}", + "deltaPeakVar^{%s}", + "n_{WT1}^{%s}", + "n_{WT2}^{%s}"}; + if (q2Bin == 7) + massConstParName = { "f_{M}^{%s}", + "#alpha_{RT1}^{%s}", + "#sigma_{RT1}^{%s}", + "#sigma_{RT2}^{%s}", + "f^{RT%s}", + "n_{RT1}^{%s}", + "#alpha_{WT1}^{%s}", + "#alpha_{WT2}^{%s}", + "#sigma_{WT1}^{%s}", + "deltaPeakVar^{%s}", + "n_{WT1}^{%s}", + "n_{WT2}^{%s}"}; + + if (q2Bin < 5) + massConstParName = { "f_{M}^{%s}", + "#alpha_{RT1}^{%s}", + "#alpha_{RT2}^{%s}", + "#sigma_{RT1}^{%s}", + "n_{RT1}^{%s}", + "n_{RT2}^{%s}", + "#alpha_{WT1}^{%s}", + "#alpha_{WT2}^{%s}", + "#sigma_{WT1}^{%s}", + "deltaPeakVar^{%s}", + "n_{WT1}^{%s}", + "n_{WT2}^{%s}"}; + + static const int nMassConstPars = massConstParName.size(); + + if (verbose) std::cout << "nMassPars = " << nMassPars << "\t nMassConstPars = " << nMassConstPars << std::endl; + + static const int nSPars = 6; + string sParName[nSPars] = {"Fs","As","A4s","A5s","A7s","A8s"}; + + double PDG_Fl [2] = {0.571,0.463}; + double PDG_Fleh[2] = {0.007,0.028}; + double PDG_Flel[2] = {0.007,0.040}; + + double bestParSim[nPars]; + double lowParSim [nPars]; + double highParSim[nPars]; + double bestPar[3][nPars]; + double lowPar [3][nPars]; + double highPar[3][nPars]; + + double bestSParSim[nSPars]; + double lowSParSim [nSPars]; + double highSParSim[nSPars]; + double bestSPar[3][nSPars]; + double lowSPar [3][nSPars]; + double highSPar[3][nSPars]; + + double bestMassParSim[3][nMassPars]; + double errMassParSim [3][nMassPars]; + double bestMassConstParSim[3][nMassConstPars]; + double errMassConstParSim [3][nMassConstPars]; + double bestMassPar[3][nMassPars]; + double errMassPar [3][nMassPars]; + double bestMassConstPar[3][nMassConstPars]; + double errMassConstPar [3][nMassConstPars]; + + double massConstVal[3][nMassConstPars]; + double massConstErr[3][nMassConstPars]; + + + + // INPUTS + // Now retrieves the angular parameters from the root file (they are saved in the fitResultsTree) + // for the full stat results + string finName = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_%s_b%i-XGBv8_unbl4.root"; + // for data subsample +// string finName = "/eos/cms/store/user/fiorendi/p5prime/angularFits/dataCR/signal-stat/simFitResult_data_fullAngularMass_Swave_%s_b5stat-49_b%i-XGBv8_unbl4.root"; + // which single year fits are available + std::vector haveSingleFits = {}; + int nYears = (int)haveSingleFits.size(); + + // Here the files used to get the constraints + string finName2 = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/massFits-LMNR-XGBv8/%s.root"; + if (q2Bin==4) finName2 = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/massFits-Jpsi-XGBv8/%s.root"; + else if (q2Bin==6) finName2 = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/massFits-Psi-XGBv8/%s.root"; + + + // Get the angular pars from the root file + auto finSim = TFile::Open(Form(finName.c_str(),"201620172018",q2Bin)); + if ( !finSim || finSim->IsZombie() ) { + cout<<"Sim file is problematic!"<Get("fitResultsTree"); +// if ( !tinSim || tinSim->IsZombie() || tinSim->GetEntries() < 1 ) { +// cout<<"Sim tree is problematic!"<Close(); +// return; +// } +// for (int iPar=0; iParSetBranchAddress((parName[iPar]+"_best").c_str(),&bestParSim[iPar]); +// tinSim->SetBranchAddress((parName[iPar]+"_low" ).c_str(),&lowParSim [iPar]); +// tinSim->SetBranchAddress((parName[iPar]+"_high").c_str(),&highParSim[iPar]); +// } +// tinSim->GetEntry(0); + + auto wsSim = (RooWorkspace*)finSim->Get("wsp_out"); + if ( !wsSim || wsSim->IsZombie() ) { + cout<<"Sim workspace is problematic!"<var(Form(massParName[iPar/3].c_str(),year[iPar%3].c_str())); + std::cout << Form(massParName[iPar/3].c_str(),year[iPar%3].c_str()) << std::endl; + bestMassParSim[iPar%3][iPar/3] = par->getValV(); + errMassParSim [iPar%3][iPar/3] = par->getError(); + if (verbose) std::cout << "sim results (free): \t" << Form(massParName[iPar/3].c_str(),year[iPar%3].c_str()) << " = " << bestMassParSim[iPar%3][iPar/3] << std::endl; + } + for (int iPar=0; iPar<3*nMassConstPars; ++iPar) { + auto par = (RooRealVar*)wsSim->var(Form(massConstParName[iPar/3].c_str(),year[iPar%3].c_str())); + std::cout << Form(massConstParName[iPar/3].c_str(),year[iPar%3].c_str()) << std::endl; + bestMassConstParSim[iPar%3][iPar/3] = par->getValV(); + errMassConstParSim [iPar%3][iPar/3] = par->getError(); + if (verbose) std::cout << "sim results (constr): \t" << Form(massConstParName[iPar/3].c_str(),year[iPar%3].c_str()) << " = " << bestMassConstParSim[iPar%3][iPar/3] << std::endl; + } + + finSim->Close(); + + + // Get pars from the single year fits, if any + for (int iYear=0; iYearIsZombie() ) { + cout<Get("fitResultsTree"); + if ( !tin || tin->IsZombie() || tin->GetEntries() < 1 ) { + cout<Close(); + return; + } + // Get P-wave angular pars + for (int iPar=0; iParSetBranchAddress((parName[iPar]+"_best").c_str(),&bestPar[iYear][iPar]); + tin->SetBranchAddress((parName[iPar]+"_low" ).c_str(),&lowPar [iYear][iPar]); + tin->SetBranchAddress((parName[iPar]+"_high").c_str(),&highPar[iYear][iPar]); + } + tin->GetEntry(0); + + // Get S-wave and mass pars + auto ws = (RooWorkspace*)fin->Get("wsp_out"); + if ( !ws || ws->IsZombie() ) { + cout<var(sParName[iPar].c_str()); + bestSPar[iYear][iPar] = par->getValV(); + lowSPar [iYear][iPar] = par->getErrorLo(); + highSPar[iYear][iPar] = par->getErrorHi(); + } + for (int iPar=0; iParvar(Form(massParName[iPar].c_str(),year[iYear].c_str())); + if (!par) { + cout<<"Error: "<getValV(); + errMassPar [iYear][iPar] = par->getError(); + if (verbose) std::cout << "single fit (free): " << Form(massParName[iPar].c_str(),year[iYear].c_str()) << " = " << bestMassPar[iYear][iPar] << std::endl; + } + for (int iPar=0; iParvar(Form(massConstParName[iPar].c_str(),year[iYear].c_str())); + if (!par) { + cout<<"Error: "<getValV(); + errMassConstPar [iYear][iPar] = par->getError(); + if (verbose) std::cout << "single fit (constr): " << Form(massConstParName[iPar].c_str(),year[iYear].c_str()) << " = " << bestMassConstPar[iYear][iPar] << std::endl; + } + + fin->Close(); + } + + // Now read constraints from the mass-pdf files + for (int iYear=0; iYear<3; ++iYear) { + auto fin2 = TFile::Open(Form(finName2.c_str(),year[iYear].c_str())); + if ( !fin2 || fin2->IsZombie() ) { + cout<Get("w"); + if ( !ws2 || ws2->IsZombie() ) { + cout<loadSnapshot(Form("reference_fit_RT_%i",q2Bin)); + double mean_rt_val = ws2->var(Form("mean_{RT}^{%i}", q2Bin))->getVal(); + double mean_rt_err = ws2->var(Form("mean_{RT}^{%i}", q2Bin))->getError(); + + for (int iPar=0; iParloadSnapshot(Form("reference_fit_WT_%i",q2Bin)); + + if ((massConstParName[iPar].find("deltaPeak") != string::npos)) { + double mean_wt_val = ws2->var(Form("mean_{WT}^{%i}", q2Bin))->getVal(); + double mean_wt_err = ws2->var(Form("mean_{WT}^{%i}", q2Bin))->getError(); + massConstVal[iYear][iPar] = mean_rt_val-mean_wt_val; + massConstErr[iYear][iPar] = sqrt(mean_rt_err*mean_rt_err+mean_wt_err*mean_wt_err); + } else if (iPar>0) { + auto par2 = (RooRealVar*)ws2->var(Form(massConstParName[iPar].c_str(),Form("%i",q2Bin))); + if ( !par2 || par2->IsZombie() ) { + cout<getValV(); + massConstErr[iYear][iPar] = par2->getError(); + if (q2Stat>=0) massConstErr[iYear][iPar] *= constr_width_fac[atoi(year[iYear].c_str())][2*q2Stat+iPar/8]; + } else { + double nrt_mc = ws2->var(Form("nRT_%i",q2Bin))->getVal(); + double nwt_mc = ws2->var(Form("nWT_%i",q2Bin))->getVal(); + massConstVal[iYear][iPar] = 1.; + double fraction = nwt_mc / (nrt_mc + nwt_mc); +// massConstErr[iYear][iPar] = fM_sigmas[atoi(year[iYear].c_str())][q2Bin] / fraction/(1-fraction); + massConstErr[iYear][iPar] = fM_sigmas[atoi(year[iYear].c_str())][q2Stat<0?q2Bin:q2Stat] / fraction/(1-fraction); + } + } + fin2->Close(); + } + + + + // Plot mass parameters + double xMass[3*(nMassPars+nMassConstPars)]; + double xMasse[3*(nMassPars+nMassConstPars)]; + + double yMassSim [3*(nMassPars+nMassConstPars)]; + double yMasseSim[3*(nMassPars+nMassConstPars)]; + double yMass [3*(nMassPars+nMassConstPars)]; + double yMasse[3*(nMassPars+nMassConstPars)]; + double yMassConst [3*(nMassPars+nMassConstPars)]; + double yMasseConst[3*(nMassPars+nMassConstPars)]; + + for (int iPar=0; iPar<3*(nMassPars+nMassConstPars); ++iPar) { + xMass[iPar] = 0.5+iPar; + xMasse[iPar] = 0.5; + yMassSim [iPar] = 0.; + yMasseSim[iPar] = 1.; + } + + for (int iPar=0; iParGetXaxis()->SetBinLabel(iBin*3 + iYear+1,Form(massParName[iBin].c_str(), this_year.c_str())); + } + } + for (int iBin=0; iBin<1; ++iBin){ + for (int iYear=0; iYear<3; ++iYear){ + string this_year = year[iYear]; + hLabM->GetXaxis()->SetBinLabel(iBin*3 + iYear+1 +3*nMassPars,Form("R^{%s}", this_year.c_str())); + } + } + for (int iBin=1; iBinGetXaxis()->SetBinLabel(iBin*3 + iYear+1 +3*nMassPars,Form(massConstParName[iBin].c_str(), this_year.c_str())); + if (massConstParName[iBin].find("eltaPeak") !=std::string::npos) + hLabM->GetXaxis()->SetBinLabel(iBin*3 + iYear+1 +3*nMassPars,Form("#Deltam^{%s}", this_year.c_str())); + } + } + double yRangeM = 4; + if (q2Bin==4) yRangeM = 30;//30; + hLabM->SetMinimum(-1*yRangeM); + hLabM->SetMaximum(yRangeM); + hLabM->Draw(); + hLabM->GetXaxis()->LabelsOption("v"); + + + auto grSimM = new TGraphErrors(3*(nMassPars+nMassConstPars),xMass,yMassSim,xMasse,yMasseSim); + grSimM->SetName("grSimM"); + grSimM->SetFillColor(30); + grSimM->Draw("2"); + + TLine lineM (0,0,3*(nMassPars+nMassConstPars),0); + lineM.SetLineWidth(2); + lineM.SetLineStyle(2); + lineM.SetLineColor(13); + lineM.Draw(); + +// auto grM = new TGraphErrors(3*(nMassPars+nMassConstPars),xMass,yMass,xMasse,yMasse); +// grM->SetName("grM"); +// grM->SetLineColor(2); +// grM->SetLineWidth(2); +// grM->Draw("P"); +// + auto grConst = new TGraphErrors(3*(nMassPars+nMassConstPars),xMass,yMassConst,xMasse,yMasseConst); + grConst->SetName("grConst"); + grConst->SetLineColor(1);\ + grConst->SetLineWidth(2); + grConst->SetFillStyle(0); + grConst->SetMarkerStyle(20); + grConst->Draw("5p"); + + TLegend legM (0.1,0.73,0.35,0.9); + legM.AddEntry(grSimM,"Simultaneous result","f"); +// legM.AddEntry(grM,"Single-year result","lep"); + legM.AddEntry(grConst,"Constraint","fp"); + legM.Draw(); + + canvM.SaveAs(("plotSimFit4d_d/comparisonSimFit_massParameters_"+shortString + statString + "_XGBv8.pdf").c_str()); + + // Plot contrained mass parameters with traditional pull plots + { + double xMass[3*nMassConstPars]; + double xMasse[3*nMassConstPars]; + + double yMassSim [3*nMassConstPars]; + double yMasseSim[3*nMassConstPars]; + double yMass [3*nMassConstPars]; + double yMasse[3*nMassConstPars]; + double yMassConst [3*nMassConstPars]; + double yMasseConst[3*nMassConstPars]; + + for (int iPar=0; iPar<3*nMassConstPars; ++iPar) { + xMass[iPar] = 0.5+iPar; + xMasse[iPar] = 0.5; + yMassConst [iPar] = 0.; + yMasseConst[iPar] = 1.; + } + for (int iPar=0; iParSetGridx(); + + string plotTitleM = Form("Pulls of constrained fit parameters - q^{2} bin %i%s;(#theta_{fit} - #theta_{constr})/#sigma(#theta)_{constr};",q2Bin, q2Stat>=0?Form(" subsample (q2 bin %i stat)",q2Stat):""); + + double yRangeM = 4; + if (q2Bin==4) yRangeM = 5; + auto hLabM = new TH2S ("hLabPull",plotTitleM.c_str(), + 1,-1*yRangeM,yRangeM, + 3*nMassConstPars,0,3*nMassConstPars); + for (int iBin=0; iBin<3; ++iBin) + hLabM->GetYaxis()->SetBinLabel(iBin+1,Form("R^{%s}",year[iBin%3].c_str())); + for (int iBin=3; iBin<3*nMassConstPars; ++iBin){ + hLabM->GetYaxis()->SetBinLabel(iBin+1,Form(massConstParName[iBin/3].c_str(),year[iBin%3].c_str())); + if (massConstParName[iBin/3].find("eltaPeak") !=std::string::npos) +// hLabM->GetYaxis()->SetBinLabel(iBin+1, "#Deltam"); + hLabM->GetYaxis()->SetBinLabel(iBin+1,Form("#Deltam^{%s}",year[iBin%3].c_str())); + } + hLabM->Draw(); + + auto grSimM = new TGraphErrors(3*nMassConstPars,yMassSim,xMass,yMasseSim,xMasse); + grSimM->SetName("grSimPull"); + grSimM->SetLineColor(30); + grSimM->SetLineWidth(2); + grSimM->Draw("P"); + + canvM.SaveAs(("plotSimFit4d_d/comparisonSimFit_massParametersPulls_"+shortString+ statString +"_XGBv8.pdf").c_str()); + + } + +} diff --git a/plotMultiDataFit.cc b/plotMultiDataFit.cc index 93c7101..dbae017 100644 --- a/plotMultiDataFit.cc +++ b/plotMultiDataFit.cc @@ -8,7 +8,7 @@ using namespace std; static const int nBins = 8; -float binBorders [nBins+1] = { 1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; +float binBorders [nBins+1] = { 1.1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; static const int nUncBins = 100; double binsUnc [nUncBins+1]; @@ -71,11 +71,10 @@ void plotMultiDataFit () vq2Bins.push_back(q2stat); TChain fitResultsTree ("fitResultsTree",""); - string filename = Form("/eos/cms/store/group/phys_bphys/fiorendi/p5prime/UML-fit/simFitResults4d/simFitResult_data_fullAngularMass_Swave_201620172018_b%istat-*_b%i-XGBv8.root",q2stat,q2Bin); + string filename = Form("/eos/cms/store/user/fiorendi/p5prime/angularFits/dataCR/signal-stat/simFitResult_data_fullAngularMass_Swave_201620172018_b%istat-*_b%i-XGBv8_unbl4.root",q2stat,q2Bin); fitResultsTree.Add(filename.c_str()); -// string filename_fR = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%ip1_XGBv8.root",q2Bin); // link will be updated - string filename_fR = Form("/afs/cern.ch/work/d/dini/public/ResonantXGBv8/Corr_frac_sigma/simFitResults4d/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8.root",q2Bin); + string filename_fR = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); TFile* filein_fR = TFile::Open(filename_fR.c_str()); TTree* fitResultsTree_fR = (TTree*)filein_fR->Get("fitResultsTree"); if (!fitResultsTree_fR || fitResultsTree_fR->GetEntries() != 1) { @@ -216,8 +215,10 @@ void plotMultiDataFit () } legUnc->SetBorderSize(0); - legUnc->AddEntry(vHistErrH[iPar][0],Form("%i higher",vq2Bins[0]),"f"); - legUnc->AddEntry(vHistErrL[iPar][0],Form("%i lower",vq2Bins[0]),"f"); + legUnc->AddEntry(vHistErrH[iPar][0],Form("%i higher, mean: %.3f",vq2Bins[0], vHistErrH[iPar][0]->GetMean()),"f"); + legUnc->AddEntry(vHistErrL[iPar][0],Form("%i lower, mean: %.3f",vq2Bins[0], vHistErrL[iPar][0]->GetMean()),"f"); + + printf("%s:\t Bin 0 Mean (high/low) = %.5f\t\n",parName[iPar].c_str(),(vHistErrH[iPar][0]->GetMean()+vHistErrL[iPar][0]->GetMean())/2); for (int iBin=1; iBincd(); @@ -239,8 +240,10 @@ void plotMultiDataFit () leg->AddEntry(vHistBest[iPar][iBin],Form("%i [Bias:%.3f RMS:%.3f]",vq2Bins[iBin],vBias[iPar][iBin],vRMS[iPar][iBin]),"l"); legPull->AddEntry(vHistPull[iPar][iBin],Form("%i [Mean:%.3f RMS:%.3f]",vq2Bins[iBin],vHistPull[iPar][iBin]->GetMean(),vHistPull[iPar][iBin]->GetRMS()),"l"); - legUnc->AddEntry(vHistErrH[iPar][iBin],Form("%i higher",vq2Bins[iBin]),"f"); - legUnc->AddEntry(vHistErrL[iPar][iBin],Form("%i lower",vq2Bins[iBin]),"f"); + legUnc->AddEntry(vHistErrH[iPar][iBin],Form("%i higher, mean: %.3f",vq2Bins[iBin], vHistErrH[iPar][iBin]->GetMean()),"f"); + legUnc->AddEntry(vHistErrL[iPar][iBin],Form("%i lower, mean: %.3f",vq2Bins[iBin], vHistErrL[iPar][iBin]->GetMean()),"f"); + + printf("%s:\tBin 5 Mean (high/low) = %.5f\t\n",parName[iPar].c_str(),(vHistErrH[iPar][iBin]->GetMean()+vHistErrL[iPar][iBin]->GetMean())/2); } vHistBest[iPar][0]->GetYaxis()->SetRangeUser(0,1.1*ymax); diff --git a/plotMultiFit.cc b/plotMultiFit.cc index e218483..51620f8 100644 --- a/plotMultiFit.cc +++ b/plotMultiFit.cc @@ -5,20 +5,14 @@ #include #include +#include "plot_utils.h" + using namespace std; -static const int nBins = 8; -float binBorders [nBins+1] = { 1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; static const int nUncBins = 100; double binsUnc [nUncBins+1]; -static const int nPars = 8; -string parName [nPars] = {"Fl","P1","P2","P3","P4p","P5p","P6p","P8p"}; -string parTitle[nPars] = {"F_{L}","P_{1}","P_{2}","P_{3}","P'_{4}","P'_{5}","P'_{6}","P'_{8}"}; -double parMin [nPars] = {0,-1,-0.5,-0.5,-1*sqrt(2),-1*sqrt(2),-1*sqrt(2),-1*sqrt(2)}; -double parMax [nPars] = {1, 1, 0.5, 0.5, sqrt(2), sqrt(2), sqrt(2), sqrt(2)}; - static const int nQuant = 4; double quantPerc [nQuant] = {0.025,0.16,0.84,0.975}; @@ -28,14 +22,15 @@ int colors [12] = { 633, 417, 879, 857, 839, 801, 921, 607, 807, 419, 907, 402 } double diffMax = 0.0999; // double diffMax = 0.0499; +float statUnc[nPars][nBins]; void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref4dFit = false) { gROOT->SetBatch(true); - // whichSamples=0 -> plot fit results to 3D MC subsamples - // whichSamples=1 -> plot fit results to 4D MC subsamples - // whichSamples=2 -> plot fit results to 4D MC subsamples + toy background + // whichSamples=0 -> plot fit results of 3D MC subsamples + // whichSamples=1 -> plot fit results of 4D MC subsamples + // whichSamples=2 -> plot fit results of 4D MC subsamples + toy background string fitName = "3D fits to data-like signal MC subsamples"; if (whichSamples==1) fitName = "4D fits to data-like signal MC subsamples"; if (whichSamples==2) fitName = "4D fits to signal MC + toy bkg samples"; @@ -56,6 +51,7 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref vector vq2Bins (0); vector table (0); + vector stat_unc (0); binsUnc[0]=0.006; for (int i=0; iGet("fitResultsTree"); if (!fitResultsTree_fR || fitResultsTree_fR->GetEntries() != 1) { @@ -112,17 +112,17 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref for (int iPar=0; iParSetBranchAddress(Form("%s_best",parName[iPar].c_str()),&vBest[iPar]); - fitResultsTree_fR->SetBranchAddress(Form("%s_high",parName[iPar].c_str()),&vHigh[iPar]); - fitResultsTree_fR->SetBranchAddress(Form("%s_low" ,parName[iPar].c_str()),&vLow [iPar]); + fitResultsTree_fR->SetBranchAddress(Form("%s_best",ParName[iPar].c_str()),&vBest[iPar]); + fitResultsTree_fR->SetBranchAddress(Form("%s_high",ParName[iPar].c_str()),&vHigh[iPar]); + fitResultsTree_fR->SetBranchAddress(Form("%s_low" ,ParName[iPar].c_str()),&vLow [iPar]); - vHistBest[iPar].push_back( new TH1D(Form("hBest%i%i",q2Bin,iPar),Form("%s results of data-like MC sample fits - q2 bin %i;%s;# of results",parTitle[iPar].c_str(),q2Bin,parTitle[iPar].c_str()),100,parMin[iPar],parMax[iPar]) ); - vHistErrH[iPar].push_back( new TH1D(Form("hErrH%i%i",q2Bin,iPar),Form("%s MINOS uncertainties of data-like MC sample fits - q2 bin %i;#sigma(%s);# of results",parTitle[iPar].c_str(),q2Bin,parTitle[iPar].c_str()),nUncBins,binsUnc) ); - vHistErrL[iPar].push_back( new TH1D(Form("hErrL%i%i",q2Bin,iPar),Form("%s MINOS uncertainties of data-like MC sample fits - q2 bin %i;#sigma(%s);# of results",parTitle[iPar].c_str(),q2Bin,parTitle[iPar].c_str()),nUncBins,binsUnc) ); + vHistBest[iPar].push_back( new TH1D(Form("hBest%i%i",q2Bin,iPar),Form("%s results of data-like MC sample fits - q2 bin %i;%s;# of results",ParTitle[iPar].c_str(),q2Bin,ParTitle[iPar].c_str()),100,ParMin[iPar],ParMax[iPar]) ); + vHistErrH[iPar].push_back( new TH1D(Form("hErrH%i%i",q2Bin,iPar),Form("%s MINOS uncertainties of data-like MC sample fits - q2 bin %i;#sigma(%s);# of results",ParTitle[iPar].c_str(),q2Bin,ParTitle[iPar].c_str()),nUncBins,binsUnc) ); + vHistErrL[iPar].push_back( new TH1D(Form("hErrL%i%i",q2Bin,iPar),Form("%s MINOS uncertainties of data-like MC sample fits - q2 bin %i;#sigma(%s);# of results",ParTitle[iPar].c_str(),q2Bin,ParTitle[iPar].c_str()),nUncBins,binsUnc) ); vHistBest[iPar].back()->SetLineColor(colors[iColor]); vHistErrH[iPar].back()->SetLineColor(colors[iColor]); @@ -159,6 +159,7 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref string firstLine = "\\textbf{$q^2$-bin}"; string line = Form("%i",q2Bin); + for (int iPar=0; iPar "<<( vRMS[iPar].back() - vBias[iPar].back() * vBias[iPar].back() ) / ( fitResultsTree.GetEntries() - 1 )<cd(); vHistBest[iPar][0]->Draw(); @@ -223,18 +229,18 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref TLegend* leg; TLegend* legUnc; - if ( parName[iPar].compare("P4p")==0 || parName[iPar].compare("P5p")==0 || parName[iPar].compare("P1")==0 ) + if ( ParName[iPar].compare("P4p")==0 || ParName[iPar].compare("P5p")==0 || ParName[iPar].compare("P1")==0 ) leg = new TLegend(0.67,0.57,0.87,0.87,"q^{2} bin"); else leg = new TLegend(0.15,0.57,0.35,0.87,"q^{2} bin"); - if ( parName[iPar].compare("P4p")==0 || parName[iPar].compare("P8p")==0 || parName[iPar].compare("P3")==0 ) + if ( ParName[iPar].compare("P4p")==0 || ParName[iPar].compare("P8p")==0 || ParName[iPar].compare("P3")==0 ) legUnc = new TLegend(0.15,0.57,0.35,0.87,"q^{2} bin"); else legUnc = new TLegend(0.67,0.57,0.87,0.87,"q^{2} bin"); legUnc->SetNColumns(2); if (nPlotBins>1) { - vHistBest[iPar][0]->SetTitle( Form("%s results of %s",parTitle[iPar].c_str(),fitName.c_str()) ); - vHistErrH[iPar][0]->SetTitle( Form("%s MINOS uncertainties of %s",parTitle[iPar].c_str(),fitName.c_str()) ); + vHistBest[iPar][0]->SetTitle( Form("%s results of %s",ParTitle[iPar].c_str(),fitName.c_str()) ); + vHistErrH[iPar][0]->SetTitle( Form("%s MINOS uncertainties of %s",ParTitle[iPar].c_str(),fitName.c_str()) ); leg->SetBorderSize(0); leg->AddEntry(vHistBest[iPar][0],Form("%i [Bias:%.3f RMS:%.3f]",vq2Bins[0],vBias[iPar][0],vRMS[iPar][0]),"l"); } @@ -290,8 +296,8 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref if (ref4dFit) toyConfString = toyConfString + "_vs4DfullMC"; else toyConfString = toyConfString + "_vs3DfullMC"; - cDistr[iPar]->SaveAs(Form("plotSimFit_d/simfit_recoMC_%s_dist-%s_p%i.pdf",toyConfString.c_str(),parName[iPar].c_str(),parity)); - cUncert[iPar]->SaveAs(Form("plotSimFit_d/simfit_recoMC_%s_uncert-%s_p%i.pdf",toyConfString.c_str(),parName[iPar].c_str(),parity)); + cDistr[iPar]->SaveAs(Form("plotSimFit_d/simfit_recoMC_%s_dist-%s_p%i.pdf",toyConfString.c_str(),ParName[iPar].c_str(),parity)); + cUncert[iPar]->SaveAs(Form("plotSimFit_d/simfit_recoMC_%s_uncert-%s_p%i.pdf",toyConfString.c_str(),ParName[iPar].c_str(),parity)); // Plot resutls vs q2 @@ -347,15 +353,14 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref aQuantInnerCenter[vq2Bins[iBin]] = 0.5 * ( quantVal[1] + quantVal[2] ); aQuantOuterCenter[vq2Bins[iBin]] = 0.5 * ( quantVal[0] + quantVal[3] ); aQuantInnerError[vq2Bins[iBin]] = 0.5 * fabs( quantVal[1] - quantVal[2] ); - std::cout << "bin " << iBin << " aQuantInnerError[vq2Bins]: " << aQuantInnerError[vq2Bins[iBin]] << std::endl; aQuantOuterError[vq2Bins[iBin]] = 0.5 * fabs( quantVal[0] - quantVal[3] ); } auto GrReco = new TGraphAsymmErrors(nBins, q2Val, aReco, q2Err, q2Err, aRecoErrL, aRecoErrH); GrReco->SetName(Form("GrReco%i",iPar)); - GrReco->SetTitle(Form("%s results from %s",parTitle[iPar].c_str(),fitName.c_str())); - GrReco->GetYaxis()->SetTitle(parTitle[iPar].c_str()); + GrReco->SetTitle(Form("%s results from %s",ParTitle[iPar].c_str(),fitName.c_str())); + GrReco->GetYaxis()->SetTitle(ParTitle[iPar].c_str()); auto Gr = new TGraphErrors(nBins, q2Val, aMean, q2Err, aMeanErr); Gr->SetName(Form("Gr%i",iPar)); auto GrDiff = new TGraphAsymmErrors(nBins, q2Val, aBias, q2Err, q2Err, aBiasErrL, aBiasErrH); @@ -383,20 +388,8 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref GrQuantOut->SetFillStyle(0); // Grey bands for resonant regions - double ResX [2] = {0.5*(binBorders[5]+binBorders[4]),0.5*(binBorders[7]+binBorders[6])}; - double ResXe[2] = {0.5*(binBorders[5]-binBorders[4]),0.5*(binBorders[7]-binBorders[6])}; - double ResY [2] = {0.5*(parMax[iPar]+parMin[iPar]),0.5*(parMax[iPar]+parMin[iPar])}; - double ResYe[2] = {0.498*(parMax[iPar]-parMin[iPar]),0.498*(parMax[iPar]-parMin[iPar])}; - double ResD [2] = {0,0}; - double ResDe[2] = { 0.98*diffMax, 0.98*diffMax}; - TGraphErrors *resCover = new TGraphErrors(2,ResX,ResY,ResXe,ResYe); - resCover->SetName(Form("resCover%i",iPar)); - resCover->SetFillColor(18); - resCover->SetFillStyle(1001); - TGraphErrors *resDiffCover = new TGraphErrors(2,ResX,ResD,ResXe,ResDe); - resDiffCover->SetName(Form("resDiffCover%i",iPar)); - resDiffCover->SetFillColor(18); - resDiffCover->SetFillStyle(1001); + TGraphErrors *resCover = getGreyResonanceGraph(iPar); + TGraphErrors *resDiffCover = getGreyResonancePull(iPar, diffMax); // Legend TLegend *legRes; @@ -429,13 +422,13 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref GrReco->GetYaxis()->SetLabelSize(0.); GrReco->GetYaxis()->SetTitleSize(0.); - GrReco->GetYaxis()->SetRangeUser(parMin[iPar],parMax[iPar]); - TGaxis *axis = new TGaxis( GrReco->GetXaxis()->GetXmin(), parMin[iPar]+0.01, - GrReco->GetXaxis()->GetXmin(), parMax[iPar], - parMin[iPar]+0.01,parMax[iPar], + GrReco->GetYaxis()->SetRangeUser(ParMin[iPar],ParMax[iPar]); + TGaxis *axis = new TGaxis( GrReco->GetXaxis()->GetXmin(), ParMin[iPar]+0.01, + GrReco->GetXaxis()->GetXmin(), ParMax[iPar], + ParMin[iPar]+0.01,ParMax[iPar], 510, ""); axis->SetName(Form("axis%i",iPar)); - axis->SetTitle(parTitle[iPar].c_str()); + axis->SetTitle(ParTitle[iPar].c_str()); axis->SetLabelFont(43); axis->SetLabelSize(20); axis->Draw(); @@ -485,7 +478,7 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref line->Draw(); resDiffCover->Draw("e2"); - cResult[iPar]->SaveAs(Form("plotSimFit_d/simfit_recoMC_%s_results-%s_p%i.pdf",toyConfString.c_str(),parName[iPar].c_str(),parity)); + cResult[iPar]->SaveAs(Form("plotSimFit_d/simfit_recoMC_%s_results-%s_p%i.pdf",toyConfString.c_str(),ParName[iPar].c_str(),parity)); } @@ -493,4 +486,12 @@ void plotMultiFit (int binIndex=-1, int parity=1, int whichSamples = 2, bool ref for (int iLine=0; iLineIsZombie() ) { diff --git a/plotSimFitComparison_manual.cc b/plotSimFitComparisonFromLogs.cc similarity index 89% rename from plotSimFitComparison_manual.cc rename to plotSimFitComparisonFromLogs.cc index 4b4039b..4d76f4d 100644 --- a/plotSimFitComparison_manual.cc +++ b/plotSimFitComparisonFromLogs.cc @@ -2,29 +2,31 @@ string year[3] = {"2016","2017","2018"}; -void plotSimFitComparison_manual(int q2Bin = 4) +void plotSimFitComparisonFromLogs(int q2Bin = 6) { gROOT->SetBatch(true); bool verbose = false; - string shortString = Form("b%ip1",q2Bin); + string shortString = Form("b%ip0",q2Bin); static const int nPars = 8; string parName[nPars] = {"Fl","P1","P2","P3","P4p","P5p","P6p","P8p"}; - static const int nMassPars = 2; - string massParName[nMassPars] = {"mean_{RT}^{%s}", - "fsig_"+shortString+"_%s"}; +// static const int nMassPars = 2; + string massParName[] = {"mean_{RT}^{%s}", + "fsig_"+shortString+"_%s"}; // "slope^{%s}"}; + static const int nMassPars = sizeof(massParName)/sizeof(string); + static const int nMassConstPars = 14; string massConstParName[nMassConstPars] = {"f_{M}^{%s}", "#alpha_{RT1}^{%s}", "#alpha_{RT2}^{%s}", "#sigma_{RT1}^{%s}", - "#sigma_{RT2}^{%s}", - "f^{RT%s}", + "#sigma_{RT2}^{%s}", // not for jpsi + "f^{RT%s}", "n_{RT1}^{%s}", "n_{RT2}^{%s}", "#alpha_{WT1}^{%s}", @@ -32,7 +34,7 @@ void plotSimFitComparison_manual(int q2Bin = 4) "#sigma_{WT1}^{%s}", "deltaPeakVar^{%s}", "n_{WT1}^{%s}", - "n_{WT2}^{%s}"}; + "n_{WT2}^{%s}"}; static const int nSPars = 6; string sParName[nSPars] = {"Fs","As","A4s","A5s","A7s","A8s"}; @@ -69,12 +71,15 @@ void plotSimFitComparison_manual(int q2Bin = 4) // First reads values of parameters which are not P_i or F_L from the log file of the simultaneous fit - string finLogName = "/afs/cern.ch/work/d/dini/public/ResonantXGBv8/Corr_frac_sigma/Bin%i-allYears-Bern.log"; - // string finLogName = "logs_parSub/simfit_data_fullAngularMass_Swave_0_%i_0.out"; - // string finLogName = "logs_simFit4d/simfit_data_fullAngularMass_Swave_%i_1_0_0_0_0_1_2016_2017_2018.log"; + // for AN v12 + // string finLogName = "/afs/cern.ch/work/d/dini/public/ResonantSwapCut/Bin4SimBern/logs_simFit4d/simfit_data_fullAngularMass_Swave_%i_0_0_0_0_1_8_0_1_2_2016_2017_2018_unbl4.log"; + // for AN v13 + string finLogName = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simfit_data_fullAngularMass_Swave_%i_0_0_0_0_1_8_0_1_2_2016_2017_2018_unbl4.log"; + if (q2Bin == 6) + finLogName = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simfit_data_fullAngularMass_Swave_%i_0_0_0_0_0_8_0_1_2_2016_2017_2018_unbl4.log"; ifstream fin_log(Form(finLogName.c_str(),q2Bin)); if ( !fin_log ) { - cout<<"Log file is problematic!"< haveSingleFits = {2016,2017,2018}; int nYears = (int)haveSingleFits.size(); @@ -141,6 +148,9 @@ void plotSimFitComparison_manual(int q2Bin = 4) // Get pars from the single year fits, if any for (int iYear=0; iYearIsZombie() ) { cout<IsZombie() ) { - cout<Get("w"); @@ -212,9 +222,12 @@ void plotSimFitComparison_manual(int q2Bin = 4) ws2->loadSnapshot(Form("reference_fit_RT_%i",q2Bin)); double mean_rt_val = ws2->var(Form("mean_{RT}^{%i}", q2Bin))->getVal(); double mean_rt_err = ws2->var(Form("mean_{RT}^{%i}", q2Bin))->getError(); + for (int iPar=0; iParloadSnapshot(Form("reference_fit_WT_%i",q2Bin)); - if (iPar==11) { + if ((massConstParName[iPar].find("WT") != string::npos)) + ws2->loadSnapshot(Form("reference_fit_WT_%i",q2Bin)); + + if ((massConstParName[iPar].find("deltaPeak") != string::npos)) { double mean_wt_val = ws2->var(Form("mean_{WT}^{%i}", q2Bin))->getVal(); double mean_wt_err = ws2->var(Form("mean_{WT}^{%i}", q2Bin))->getError(); massConstVal[iYear][iPar] = mean_rt_val-mean_wt_val; @@ -223,6 +236,8 @@ void plotSimFitComparison_manual(int q2Bin = 4) auto par2 = (RooRealVar*)ws2->var(Form(massConstParName[iPar].c_str(),Form("%i",q2Bin))); if ( !par2 || par2->IsZombie() ) { cout<Print(); + cout << "from file: " << Form(finName2.c_str(),year[iYear].c_str()) << endl; return; } massConstVal[iYear][iPar] = par2->getValV(); @@ -329,7 +344,7 @@ void plotSimFitComparison_manual(int q2Bin = 4) leg.AddEntry(gr[iYear], Form("Result %i", haveSingleFits[iYear]),"lep"); leg.Draw(); - canv.SaveAs(("plotSimFit4d_d/comparisonSimFit_result_"+shortString+"_XGBv8_newMaster.pdf").c_str()); + canv.SaveAs(("plotSimFit4d_d/comparisonSimFit_result_"+shortString+"_XGBv8.pdf").c_str()); // Plot S-wave parameters double xS[nYears][nSPars]; @@ -458,6 +473,8 @@ void plotSimFitComparison_manual(int q2Bin = 4) for (int iBin=nYears; iBinGetXaxis()->SetBinLabel(iBin+1+nYears*nMassPars,Form(massConstParName[iBin/nYears].c_str(), this_year.c_str())); + if (massConstParName[iBin/nYears].find("eltaPeak") !=std::string::npos) + hLabM->GetXaxis()->SetBinLabel(iBin+1+nYears*nMassPars,Form("#Deltam^{%s}", this_year.c_str())); } double yRangeM = 9; @@ -531,7 +548,7 @@ void plotSimFitComparison_manual(int q2Bin = 4) auto vp = canvM.cd(); vp->SetGridx(); - string plotTitleM = Form("Pulls of contrained fit parameters - q2 bin %i;(#theta_{fit} - #theta_{constr})/#sigma(#theta)_{constr};",q2Bin); + string plotTitleM = Form("Pulls of constrained fit parameters - q2 bin %i;(#theta_{fit} - #theta_{constr})/#sigma(#theta)_{constr};",q2Bin); double yRangeM = 9; if (q2Bin==4) yRangeM = 20; auto hLabM = new TH2S ("hLabPull",plotTitleM.c_str(), @@ -539,8 +556,11 @@ void plotSimFitComparison_manual(int q2Bin = 4) 3*nMassConstPars,0,3*nMassConstPars); for (int iBin=0; iBin<3; ++iBin) hLabM->GetYaxis()->SetBinLabel(iBin+1,Form("R^{%s}",year[iBin%3].c_str())); - for (int iBin=3; iBin<3*nMassConstPars; ++iBin) + for (int iBin=3; iBin<3*nMassConstPars; ++iBin){ hLabM->GetYaxis()->SetBinLabel(iBin+1,Form(massConstParName[iBin/3].c_str(),year[iBin%3].c_str())); + if (massConstParName[iBin/3].find("eltaPeak") !=std::string::npos) + hLabM->GetYaxis()->SetBinLabel(iBin+1,Form("#Deltam^{%s}",year[iBin%3].c_str())); + } hLabM->Draw(); auto grSimM = new TGraphErrors(3*nMassConstPars,yMassSim,xMass,yMasseSim,xMasse); diff --git a/plotSimFitResults.cc b/plotSimFitResults.cc index 8136b14..cf89f34 100644 --- a/plotSimFitResults.cc +++ b/plotSimFitResults.cc @@ -1,6 +1,8 @@ using namespace RooFit; using namespace std; +#include "plot_utils.h" + /* code to plot results from the Simultaneous fit of multiple years and compare them to - GEN level results @@ -8,75 +10,18 @@ code to plot results from the Simultaneous fit of multiple years and compare the (should be provided) */ - -static const int nBins = 8; // was 8 -float binBorders [nBins+1] = { 1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; -// static const int nBins = 9; -// float binBorders [nBins+1] = { 1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16, 19}; - -static const int nPars = 8; -string ParName [nPars] = { "Fl", "P1", "P2", "P3", "P4p", "P5p", "P6p", "P8p" }; -string ParTitle [nPars] = { "F_{L}", "P_{1}", "P_{2}", "P_{3}", "P'_{4}", "P'_{5}", "P'_{6}", "P'_{8}" }; -double ParMin [nPars] = {0,-1,-0.5,-0.5,-sqrt(2),-sqrt(2),-sqrt(2),-sqrt(2)}; -double ParMax [nPars] = {1,1,0.5,0.5,sqrt(2),sqrt(2),sqrt(2),sqrt(2)}; float diffMax = 0.0999; - +double syst_matrix[nBins][nPars]; EColor colorlist[] = {kBlue, kRed, kOrange, kGreen, kMagenta}; - TCanvas* c[nPars]; - -void fillVectors(RooFitResult* fitResult, int ParIndx, double Res[][nBins], double ErrH[][nBins], double ErrL[][nBins], - double Diff[][nBins], double DErrH[][nBins], double DErrL[][nBins], - double genRes[nBins], double genErrH[nBins], double genErrL[nBins], - int year, int ibin){ - RooRealVar* Par = (RooRealVar*)fitResult->floatParsFinal().find(ParName[ParIndx].c_str()); - Res[year][ibin] = Par->getValV(); - ErrH[year][ibin] = Par->getErrorHi(); - ErrL[year][ibin] = -1*Par->getErrorLo(); - Diff[year][ibin] = Res[year][ibin]-genRes[ibin]; - DErrH[year][ibin] = sqrt(ErrH[year][ibin]*ErrH[year][ibin]+genErrL[ibin]*genErrL[ibin]); - DErrL[year][ibin] = sqrt(ErrL[year][ibin]*ErrL[year][ibin]+genErrH[ibin]*genErrH[ibin]); -} - -void fillVectors(RooFitResult* fitResult, int ParIndx, double Res[nBins], double ErrH[nBins], double ErrL[nBins], int ibin){ - RooRealVar* Par = (RooRealVar*)fitResult->floatParsFinal().find(ParName[ParIndx].c_str()); - Res[ibin] = Par->getValV(); - ErrH[ibin] = Par->getErrorHi(); - ErrL[ibin] = -1*Par->getErrorLo(); -} - -void setGraphProperties(TGraphAsymmErrors* Gr, int icolor, string year, int ParIndx, string type){ - - int addc = 0; - if (type.find("WT") != std::string::npos) addc = 1; - if (type.find("CT") != std::string::npos) addc = 2; - - Gr -> SetLineColor(colorlist[icolor] + addc); - Gr -> SetMarkerColor(colorlist[icolor] + addc); - Gr -> SetName(Form("Gr%s%i_%s" , type.c_str(), ParIndx, year.c_str()) ); - Gr -> SetLineWidth(2); -} - -void setAxisProperties(TH1F* auxE2){ - auxE2->SetStats(kFALSE); - auxE2->SetLineColor(1); - auxE2->GetXaxis()->SetTitle("q^{2} (GeV^{2})"); - auxE2->GetXaxis()->SetTitleSize(0.12); - auxE2->GetXaxis()->SetTitleOffset(0.95); - auxE2->GetXaxis()->SetLabelSize( 0.10); - auxE2->GetXaxis()->SetTickLength(0.1); - auxE2->GetYaxis()->SetTitle("(RECO - GEN)"); - auxE2->GetYaxis()->SetTitleSize(0.10); - auxE2->GetYaxis()->SetTitleOffset(0.45); - auxE2->GetYaxis()->SetRangeUser(-1*diffMax,diffMax); - auxE2->GetYaxis()->SetLabelSize(0.09); - auxE2->GetYaxis()->SetNdivisions(505); -} - void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool plotRECO, std::vector years, bool datalike) { + gROOT->SetBatch(true); + + string firstLine = "\\textbf{$q^2$-bin}"; + double q2Val [nBins]; double q2ErrH [nBins]; double q2ErrL [nBins]; @@ -90,7 +35,7 @@ void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool p double genRes [nBins]; double genErrH [nBins]; double genErrL [nBins]; - + // n_years = number of datasets considered + 1*simultaneous dataset string combDataString = ""; if (years.size() > 1){ @@ -140,42 +85,33 @@ void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool p std::vector GrCT, GrWT, Gr, GrDiffCT, GrDiffWT, GrDiff; -// string finGenName = "simFitResults/fitResult_genMC_penalty.root"; -// string finGenName = "../eff-KDE/fitResults/newphi/fitResult_genMC_forAN.root"; -// TFile* finGen = TFile::Open(finGenName.c_str()); -// if ( !finGen || finGen->IsZombie() ) { -// cout<<"Missing gen file: "<Get(Form("ws_b%ip%i_subs0_pow1.0",i,parity)); -// if (wspRes && !wspRes->IsZombie()) { -// RooFitResult* fitResultGen = (RooFitResult*)wspRes->obj("fitResult"); -// if (fitResultGen && !fitResultGen->IsZombie() && fitResultGen->status()==0 && fitResultGen->covQual()==3) { -// fillVectors(fitResultGen, ParIndx, genRes, genErrH, genErrL, i); -// } -// } -// } - - - TFile* finGen = TFile::Open("../eff-KDE/fitResults/newphi/fitResult_genMC_forAN.root"); + TFile* finGen = TFile::Open("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/GENFitResults/fitResult_genMC_penalty.root"); if ( !finGen || finGen->IsZombie() ) { - cout<<"Missing gen file: fitResults/fitResult_genMC.root"<Get(Form("fitResult_b%ip%i",i,parity)); - if (fitResultGen && !fitResultGen->IsZombie() && fitResultGen->status()==0 && fitResultGen->covQual()==3) { - fillVectors(fitResultGen, ParIndx, genRes, genErrH, genErrL, i); + RooWorkspace* wspRes = (RooWorkspace*)finGen->Get(Form("ws_b%ip1_s0_pow1.0",i)); + if (wspRes && !wspRes->IsZombie()) { + RooFitResult* fitResultGen = (RooFitResult*)wspRes->obj("fitResult"); + if (fitResultGen && !fitResultGen->IsZombie() && fitResultGen->status()==0 && fitResultGen->covQual()==3) { + fillVectors(fitResultGen, ParIndx, genRes, genErrH, genErrL, i); + } } } + // first fill for gen results // case for fitResults +// for (int i=0; iGet(Form("fitResult_b%ip%i",i,parity)); +// if (fitResultGen && !fitResultGen->IsZombie() && fitResultGen->status()==0 && fitResultGen->covQual()==3) { +// fillVectors(fitResultGen, ParIndx, genRes, genErrH, genErrL, i); +// } +// } + finGen->Close(); - cout << genRes[0] << endl; + cout<<"Closed gen file " <IsZombie() ) { RooFitResult* fitResultCT = (RooFitResult*)finReco->Get(Form("simFitResult_b%ip%it1",i,parity)); @@ -214,14 +151,19 @@ void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool p } } - // fill for full angular fit + // fill for full angular fit using RooFitResults +// if ( plotRECO && finFullReco && !finFullReco->IsZombie() ) { +// RooFitResult* fitResult = (RooFitResult*)finFullReco->Get(Form("simFitResult_b%ip%isubs0",i,parity)); +// if (fitResult && !fitResult->IsZombie() && fitResult->status()==0 && fitResult->covQual()==3) { +// fillVectors(fitResult, ParIndx, Res, ErrH, ErrL, Diff, DiffErrH, DiffErrL, genRes, genErrH, genErrL, iy, i); +// } +// } + // fill for full angular fit using TTree if ( plotRECO && finFullReco && !finFullReco->IsZombie() ) { - RooFitResult* fitResult = (RooFitResult*)finFullReco->Get(Form("simFitResult_b%ip%isubs0",i,parity)); - if (fitResult && !fitResult->IsZombie() && fitResult->status()==0 && fitResult->covQual()==3) { - fillVectors(fitResult, ParIndx, Res, ErrH, ErrL, Diff, DiffErrH, DiffErrL, genRes, genErrH, genErrL, iy, i); - } + TTree* ttree = (TTree*)finFullReco->Get("fitResultsTree"); + std::cout << "reading RECO results from : " << finFullReco->GetName() << std::endl; + fillVectors(ttree, ParIndx, Res, ErrH, ErrL, Diff, DiffErrH, DiffErrL, genRes, genErrH, genErrL, iy, i, -1); } - if ( plotRECO && finFullReco ) finFullReco->Close(); } @@ -231,23 +173,26 @@ void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool p // create TGraphs with fit results GrCT.push_back( new TGraphAsymmErrors(nBins, q2Val, ctRes[iy], q2ErrH, q2ErrL, ctErrL[iy], ctErrH[iy]) ); - setGraphProperties(GrCT[iy], iy, year, ParIndx, "CT"); + setGraphProperties(GrCT[iy], iy, colorlist[iy], year, ParIndx, "CT"); + +// void setGraphProperties(TGraphAsymmErrors* Gr, int icolor, int thiscolor, string label, int ParIndx, string type, int marker){ +// setGraphProperties(GrCT[ifile], ifile, colorlist[ifile], labels[ifile], ParIndx, "CT", 20); GrWT.push_back( new TGraphAsymmErrors(nBins, q2Val, wtRes[iy], q2ErrH, q2ErrL, wtErrL[iy], wtErrH[iy]) ); - setGraphProperties(GrWT[iy], iy, year, ParIndx, "WT"); +// setGraphProperties(GrWT[iy], iy, year, ParIndx, "WT"); Gr .push_back( new TGraphAsymmErrors(nBins, q2Val, Res[iy], q2ErrH, q2ErrL, ErrL[iy], ErrH[iy]) ); - setGraphProperties(Gr[iy], iy, year, ParIndx, "Full"); + setGraphProperties(Gr[iy], iy, colorlist[iy],year, ParIndx, "Full"); // create TGraphs with difference to GEN Results GrDiffCT.push_back( new TGraphAsymmErrors(nBins, q2ValReco[iy], ctDiff[iy], q2ErrHReco, q2ErrLReco, ctDiffErrL[iy], ctDiffErrH[iy]) ); - setGraphProperties(GrDiffCT[iy], iy, year, ParIndx, "DiffCT"); + setGraphProperties(GrDiffCT[iy], iy, colorlist[iy],year, ParIndx, "DiffCT"); GrDiffWT.push_back( new TGraphAsymmErrors(nBins, q2ValReco[iy], wtDiff[iy], q2ErrHReco, q2ErrLReco, wtDiffErrL[iy], wtDiffErrH[iy]) ); - setGraphProperties(GrDiffWT[iy], iy, year, ParIndx, "DiffWT"); + setGraphProperties(GrDiffWT[iy], iy, colorlist[iy],year, ParIndx, "DiffWT"); GrDiff.push_back( new TGraphAsymmErrors(nBins, q2ValReco[iy], Diff[iy], q2ErrHReco, q2ErrLReco, DiffErrL[iy], DiffErrH[iy]) ); - setGraphProperties(GrDiff[iy], iy, year, ParIndx, "Diff"); + setGraphProperties(GrDiff[iy], iy, colorlist[iy], year, ParIndx, "Diff"); } @@ -346,7 +291,7 @@ void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool p // first create axis TH1F* auxE2 = new TH1F(Form("auxE2%i",ParIndx), "", nBins, GrGen->GetXaxis()->GetXmin(), GrGen->GetXaxis()->GetXmax()); - setAxisProperties(auxE2); + setAxisProperties(auxE2, diffMax); auxE2->Draw(); // Bin lines @@ -370,11 +315,16 @@ void plotFitResultsBin(int parity, int ParIndx, bool plotCT, bool plotWT, bool p if (plotCT) confString = confString + "ctRes_"; if (plotWT) confString = confString + "wtRes_"; if (plotRECO) confString = confString + "recoRes_"; - c[ParIndx]->SaveAs( (confString + ParName[ParIndx] + "_" + years.back().c_str() + stat + ".pdf").c_str() ); + if (parity == 1) + c[ParIndx]->SaveAs( (confString + ParName[ParIndx] + "_" + years.back().c_str() + stat + ".pdf").c_str() ); + // now print out table for syst only for simultaneous fit + for (int ibin = 0; ibin < nBins; ibin++) { + syst_matrix[ibin][ParIndx] = Diff[3][ibin]; + } } -void plotSimFitResults(int parity, int ParIndx = -1, bool plotCT = true, bool plotWT = false, bool plotRECO = false, bool datalike = false) +void plotSimFitResults(int parity, int ParIndx = -1, bool plotCT = false, bool plotWT = false, bool plotRECO = true, bool datalike = false) { if ( parity<0 || parity>1 ) return; @@ -396,4 +346,28 @@ void plotSimFitResults(int parity, int ParIndx = -1, bool plotCT = true, bool pl else plotFitResultsBin(parity, ParIndx, plotCT, plotWT, plotRECO, years, datalike); + vector table (0); + string firstLine = "\\textbf{$q^2$-bin}"; + for (int ibin = 0; ibin < nBins; ibin++) { + string line = Form("%i",ibin); + for (int iPar = 0; iPar < nPars; iPar++) { + firstLine = firstLine + " & \\textbf{$" + ParTitle[iPar] + "$}"; + if (fabs(syst_matrix[ibin][iPar])<0.0005) + line = line + " & $<$ 0.001"; + else + line = line + Form(" & $\\pm %.3f$",fabs(syst_matrix[ibin][iPar])); + } + if (table.size()==0) { + firstLine = firstLine + " \\\\ \\hline"; + table.push_back(firstLine); + } + line = line + " \\\\"; + table.push_back(line); + line = line + " \\\\"; + } + + cout<<"===== Formatted bias table ========"< +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "RooDoubleCBFast.h" +#include "RooExponential.h" + +#include "utils.h" +#include "PdfSigRTMass.h" +#include "PdfSigWTMass.h" +#include "PdfSigAngMass.h" +#include "ShapeSigAng.h" + +using namespace RooFit; +using namespace std; + +static const int nBins = 9; + +TCanvas* c [4*nBins]; + +double power = 1.0; + +// static const int nRanges = 5; +// double rangeVal[nRanges+1] = {5.00,5.10,5.20,5.35,5.45,5.60}; +// string rangeName[nRanges] = {"lsb","ltail","peak","rtail","rsb"}; + +static const int nRanges = 3; +double rangeVal[nRanges+1] = {-1, -0.6, 0.2, 1.}; +string rangeName[nRanges] = {"lowCTK","mediumCTK","highCTK"}; + +void plot_simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, uint nSample, std::vector years) +{ + + RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ; + + // to speed up code + RooAbsReal::defaultIntegratorConfig()->getConfigSection("RooIntegrator1D").setRealValue("fixSteps",4) ; + + map< string, vector* > rangeComp; + for (int i=0; i {i}; + + string shortString = Form("b%ip%i",q2Bin,parity); + cout<<"Conf: "<0) stat = Form("_sub-%i",iSample); + string filename_nominal = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); + TFile* fin = new TFile(filename_nominal.c_str(),"READ"); +// simFitResults4d/simFitResult_data_fullAngularMass_Swave_201620172018_b7-XGBv8_wp90_unbl4.root + auto wsp_out = (RooWorkspace*)fin->Get("wsp_out"); + auto combData = (RooDataSet*)wsp_out->data("allcombData"); + auto simPdf = (RooSimultaneous*)wsp_out->pdf("simPdf"); + auto mass = (RooRealVar*)wsp_out->var("mass"); + auto ctK = (RooRealVar*)wsp_out->var("ctK"); + auto ctL = (RooRealVar*)wsp_out->var("ctL"); + auto phi = (RooRealVar*)wsp_out->var("phi"); + + string plotString = shortString + "_" + all_years; + if (nSample>0) plotString = plotString + Form("_s%i",nSample); + + int confIndex = 2*nBins*parity + q2Bin; + string longString = "Fit to reconstructed events"; + longString = longString + Form(parity==1?" (q2-bin %i even)":" (q2-bin %i odd) ",q2Bin); + + + // plot fit projections + c[confIndex] = new TCanvas (("c_"+shortString).c_str(),("Fit to RECO-level MC - "+longString).c_str(),2000,1500); + if (years.size()>1) c[confIndex]->Divide(4, years.size()); + else c[confIndex]->Divide(2,2); + + // now, plot fit projections in different defined ctK ranges + + // assign ranges to the ctK variable + for (int i=0; isetRange(rangeName[i].c_str(),rangeVal[i],rangeVal[i+1]); + + + for (const auto& [name, comp] : rangeComp) { + auto canv = new TCanvas (Form("canv%s",name.c_str()), + Form("canv%s",name.c_str()), + 1500,500*years.size()); + canv->Divide(3, years.size()); + + string ctKCut = Form("(ctK>%.2f && ctK<%.2f)", + rangeVal[comp->at(0)], + rangeVal[comp->at(0)+1]); + string rangeList = rangeName[comp->at(0)]; + + // cut on mass range + for (uint iCut=1; iCutsize(); ++iCut) { + ctKCut = ctKCut + Form("||(ctK>%.2f && ctK<%.2f)", + rangeVal[comp->at(iCut)], + rangeVal[comp->at(iCut)+1]); + rangeList = rangeList + "," + rangeName[comp->at(iCut)]; + } + + // loop on years + for (unsigned int iy = 0; iy < years.size(); iy++) { + year.clear(); year.assign(Form("%i",years[iy])); + + std::string frametitle = ""; + + std::vector frames; + frametitle = Form("mass projection of %s dataset, ctK range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( mass ->frame(Title(frametitle.c_str())) )); + frametitle = Form("cos#theta_{L} projection of %s dataset, ctK range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( ctL ->frame(Title(frametitle.c_str())) )); + frametitle = Form("#phi projection of %s dataset, ctK range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( phi ->frame(Title(frametitle.c_str())) )); + TLegend* leg = new TLegend (0.45,0.7,0.9,0.9); + + auto singleYearPdf = simPdf->getPdf(("data"+year+Form("_subs%d",iSample)).c_str()); + auto singleYearData = combData->reduce(("sample==sample::data"+year+Form("_subs%d",iSample)).c_str()); + + TLatex text_chi2; + text_chi2.SetTextSize(0.04); + + for (unsigned int fr = 0; fr < frames.size(); fr++){ + + singleYearData->plotOn(frames[fr], + MarkerColor(kRed+1), + LineColor(kRed+1), + Binning(40), + Cut(ctKCut.c_str()), + // Cut(("("+massCut+")&&sample==sample::data"+year+Form("_subs%d",iSample)).c_str()), + Name(("plData"+name+year).c_str())); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + Name(("plPDF"+name+year).c_str()), + NumCPU(8), + RooFit::Precision(1) + ); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(8), + Name(("plPDFbkg"+name+year).c_str()), + Components( ("bkg_pdf_"+year).c_str() ), + NumCPU(8), + RooFit::Precision(1) + ); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(880), + Name(("plPDFsig"+name+year).c_str()), + Components( ("PDF_sig_ang_mass_"+shortString+"_"+year).c_str() ), + NumCPU(8), + RooFit::Precision(1) + ); + + if (fr == 0) { + leg->AddEntry(frames[fr]->findObject(("plData"+name+year).c_str()), + "Data", "lep"); + leg->AddEntry(frames[fr]->findObject(("plPDF"+name+year ).c_str()), + "Total PDF","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFsig"+name+year ).c_str()), + "Signal","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFbkg"+name+year ).c_str()), + "Background","l"); + } + canv->cd(iy*3+fr+1); + gPad->SetLeftMargin(0.16); + frames[fr]->Draw(); + double chiSq = frames[fr]->chiSquare(("plPDF"+name+year).c_str(),("plData"+name+year).c_str()); + text_chi2.DrawLatexNDC(0.22,0.86,Form("#chi^{2}/N_{bins} = %.2f/ 40",chiSq*40)); + if (fr == 0) leg->Draw("same"); + } + } + + canv->SaveAs( ("plotSimFit4d_d/slices/simFitResult_data_" + plotString + "_" + name + "_ctkSlices.pdf").c_str() ); + + } + + return; + +} + + + +void plot_simfit_data_fullAngularMass_SwaveBin1(int q2Bin, int parity, uint nSample, std::vector years) +{ + if ( parity==-1 ) + for (parity=0; parity<2; ++parity) + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); + else + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); +} + +int main(int argc, char** argv) +{ + // q2-bin format: [0-8] for one bin + // [-1] for each bin recursively + // parity format: [0] even efficiency + // [1] odd efficiency + // [-1] for each parity recursively + + int q2Bin = -1; + int parity = -1; + + if ( argc > 1 ) q2Bin = atoi(argv[1]); + if ( argc > 2 ) parity = atoi(argv[2]); + + uint nSample = 0; + if ( argc > 3 ) nSample = atoi(argv[3]); + + std::vector years; + if ( argc > 4 && atoi(argv[4]) != 0 ) years.push_back(atoi(argv[4])); + else { + cout << "No specific years selected, using default: 2016" << endl; + years.push_back(2016); + } + if ( argc > 5 && atoi(argv[5]) != 0 ) years.push_back(atoi(argv[5])); + if ( argc > 6 && atoi(argv[6]) != 0 ) years.push_back(atoi(argv[6])); + + cout << "q2Bin " << q2Bin << endl; + cout << "parity " << parity << endl; + cout << "nSample " << nSample << endl; + cout << "years[0] " << years[0] << endl; +// cout << years[1] << endl; +// cout << years[2] << endl; + + + if ( q2Bin < -1 || q2Bin >= nBins ) return 1; + if ( parity < -1 || parity > 1 ) return 1; + + if ( q2Bin==-1 ) cout << "Running all the q2 bins" << endl; + if ( parity==-1 ) cout << "Running both the parity datasets" << endl; + + if ( q2Bin==-1 ) + for (q2Bin=0; q2Bin +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "RooDoubleCBFast.h" +#include "RooExponential.h" + +#include "utils.h" +#include "PdfSigRTMass.h" +#include "PdfSigWTMass.h" +#include "PdfSigAngMass.h" +#include "ShapeSigAng.h" + +using namespace RooFit; +using namespace std; + +static const int nBins = 9; + +TCanvas* c [4*nBins]; + +double power = 1.0; + +static const int nRanges = 3; +double rangeVal[nRanges+1] = {-1, -0.4, 0.4, 1.}; +string rangeName[nRanges] = {"lowCTL","mediumCTL","highCTL"}; + +void plot_simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, uint nSample, std::vector years) +{ + + RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ; + + // to speed up code + RooAbsReal::defaultIntegratorConfig()->getConfigSection("RooIntegrator1D").setRealValue("fixSteps",3) ; + + map< string, vector* > rangeComp; + for (int i=0; i {i}; + + string shortString = Form("b%ip%i",q2Bin,parity); + cout<<"Conf: "<0) stat = Form("_sub-%i",iSample); + string filename_nominal = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); + TFile* fin = new TFile(filename_nominal.c_str(),"READ"); + auto wsp_out = (RooWorkspace*)fin->Get("wsp_out"); + auto combData = (RooDataSet*)wsp_out->data("allcombData"); + auto simPdf = (RooSimultaneous*)wsp_out->pdf("simPdf"); + auto mass = (RooRealVar*)wsp_out->var("mass"); + auto ctK = (RooRealVar*)wsp_out->var("ctK"); + auto ctL = (RooRealVar*)wsp_out->var("ctL"); + auto phi = (RooRealVar*)wsp_out->var("phi"); + + string plotString = shortString + "_" + all_years; + if (nSample>0) plotString = plotString + Form("_s%i",nSample); + + // plot fit projections in different defined ctL ranges + // assign ranges to the ctL variable + for (int i=0; isetRange(rangeName[i].c_str(),rangeVal[i],rangeVal[i+1]); + + + for (const auto& [name, comp] : rangeComp) { + auto canv = new TCanvas (Form("canv%s",name.c_str()), + Form("canv%s",name.c_str()), + 1500,500*years.size()); + canv->Divide(3, years.size()); + + string ctLCut = Form("(ctL>%.2f && ctL<%.2f)", + rangeVal[comp->at(0)], + rangeVal[comp->at(0)+1]); + string rangeList = rangeName[comp->at(0)]; + + // cut on mass range + for (uint iCut=1; iCutsize(); ++iCut) { + ctLCut = ctLCut + Form("||(ctL>%.2f && ctL<%.2f)", + rangeVal[comp->at(iCut)], + rangeVal[comp->at(iCut)+1]); + rangeList = rangeList + "," + rangeName[comp->at(iCut)]; + } + + // loop on years + for (unsigned int iy = 0; iy < years.size(); iy++) { + year.clear(); year.assign(Form("%i",years[iy])); + + std::string frametitle = ""; + + std::vector frames; + frametitle = Form("mass projection of %s data, ctL range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( mass ->frame(Title(frametitle.c_str())) )); + frametitle = Form("cos#theta_{K} projection of %s data, ctL range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( ctK ->frame(Title(frametitle.c_str())) )); + frametitle = Form("#phi projection of %s data, ctL range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( phi ->frame(Title(frametitle.c_str())) )); + TLegend* leg = new TLegend (0.45,0.7,0.9,0.9); + + auto singleYearPdf = simPdf->getPdf(("data"+year+Form("_subs%d",iSample)).c_str()); + auto singleYearData = combData->reduce(("sample==sample::data"+year+Form("_subs%d",iSample)).c_str()); + + TLatex text_chi2; + text_chi2.SetTextSize(0.04); + + for (unsigned int fr = 0; fr < frames.size(); fr++){ + + singleYearData->plotOn(frames[fr], + MarkerColor(kRed+1), + LineColor(kRed+1), + Binning(40), + Cut(ctLCut.c_str()), + // Cut(("("+massCut+")&&sample==sample::data"+year+Form("_subs%d",iSample)).c_str()), + Name(("plData"+name+year).c_str())); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + Name(("plPDF"+name+year).c_str()), + NumCPU(12), + RooFit::Precision(1) + ); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(8), + Name(("plPDFbkg"+name+year).c_str()), + Components( ("bkg_pdf_"+year).c_str() ), + NumCPU(12), + RooFit::Precision(1) + ); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(880), + Name(("plPDFsig"+name+year).c_str()), + Components( ("PDF_sig_ang_mass_"+shortString+"_"+year).c_str() ), + NumCPU(12), + RooFit::Precision(1) + ); + + if (fr == 0) { + leg->AddEntry(frames[fr]->findObject(("plData"+name+year).c_str()), + "Data", "lep"); + leg->AddEntry(frames[fr]->findObject(("plPDF"+name+year ).c_str()), + "Total PDF","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFsig"+name+year ).c_str()), + "Signal","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFbkg"+name+year ).c_str()), + "Background","l"); + } + canv->cd(iy*3+fr+1); + gPad->SetLeftMargin(0.16); + frames[fr]->Draw(); + double chiSq = frames[fr]->chiSquare(("plPDF"+name+year).c_str(),("plData"+name+year).c_str()); + text_chi2.DrawLatexNDC(0.22,0.86,Form("#chi^{2}/N_{bins} = %.2f/ 40",chiSq*40)); + if (fr == 0) leg->Draw("same"); + } + } + + canv->SaveAs( ("plotSimFit4d_d/slices/simFitResult_data_" + plotString + "_" + name + "_ctlSlices.pdf").c_str() ); + + } + + return; + +} + + + +void plot_simfit_data_fullAngularMass_SwaveBin1(int q2Bin, int parity, uint nSample, std::vector years) +{ + if ( parity==-1 ) + for (parity=0; parity<2; ++parity) + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); + else + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); +} + +int main(int argc, char** argv) +{ + // q2-bin format: [0-8] for one bin + // [-1] for each bin recursively + // parity format: [0] even efficiency + // [1] odd efficiency + // [-1] for each parity recursively + + int q2Bin = -1; + int parity = -1; + + if ( argc > 1 ) q2Bin = atoi(argv[1]); + if ( argc > 2 ) parity = atoi(argv[2]); + + uint nSample = 0; + if ( argc > 3 ) nSample = atoi(argv[3]); + + std::vector years; + if ( argc > 4 && atoi(argv[4]) != 0 ) years.push_back(atoi(argv[4])); + else { + cout << "No specific years selected, using default: 2016" << endl; + years.push_back(2016); + } + if ( argc > 5 && atoi(argv[5]) != 0 ) years.push_back(atoi(argv[5])); + if ( argc > 6 && atoi(argv[6]) != 0 ) years.push_back(atoi(argv[6])); + + cout << "q2Bin " << q2Bin << endl; + cout << "parity " << parity << endl; + cout << "nSample " << nSample << endl; + cout << "years[0] " << years[0] << endl; +// cout << years[1] << endl; +// cout << years[2] << endl; + + + if ( q2Bin < -1 || q2Bin >= nBins ) return 1; + if ( parity < -1 || parity > 1 ) return 1; + + if ( q2Bin==-1 ) cout << "Running all the q2 bins" << endl; + if ( parity==-1 ) cout << "Running both the parity datasets" << endl; + + if ( q2Bin==-1 ) + for (q2Bin=0; q2Bin +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "RooDoubleCBFast.h" +#include "RooExponential.h" + +#include "utils.h" +#include "PdfSigRTMass.h" +#include "PdfSigWTMass.h" +#include "PdfSigAngMass.h" +#include "ShapeSigAng.h" + +using namespace RooFit; +using namespace std; + +static const int nBins = 9; + +TCanvas* c [4*nBins]; + +double power = 1.0; + +static const int nRanges = 3; +double rangeVal[nRanges+1] = {5.00,5.18,5.40,5.60}; +string rangeName[nRanges] = {"lowSB","signal","highSB"}; + +void plot_simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, uint nSample, std::vector years) +{ + + RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ; + + // to speed up code + RooAbsReal::defaultIntegratorConfig()->getConfigSection("RooIntegrator1D").setRealValue("fixSteps",4) ; + + map< string, vector* > rangeComp; + for (int i=0; i {i}; + + string shortString = Form("b%ip%i",q2Bin,parity); + cout<<"Conf: "<0) stat = Form("_sub-%i",iSample); + string filename_nominal = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); + TFile* fin = new TFile(filename_nominal.c_str(),"READ"); + auto wsp_out = (RooWorkspace*)fin->Get("wsp_out"); + auto combData = (RooDataSet*)wsp_out->data("allcombData"); + auto simPdf = (RooSimultaneous*)wsp_out->pdf("simPdf"); + auto mass = (RooRealVar*)wsp_out->var("mass"); + auto ctK = (RooRealVar*)wsp_out->var("ctK"); + auto ctL = (RooRealVar*)wsp_out->var("ctL"); + auto phi = (RooRealVar*)wsp_out->var("phi"); + + string plotString = shortString + "_" + all_years; + if (nSample>0) plotString = plotString + Form("_s%i",nSample); + + + // now, plot fit projections in different defined mass ranges + + // assign ranges to the mass variable + for (int i=0; isetRange(rangeName[i].c_str(),rangeVal[i],rangeVal[i+1]); + + + for (const auto& [name, comp] : rangeComp) { + auto canv = new TCanvas (Form("canv%s",name.c_str()), + Form("canv%s",name.c_str()), + 1500,500*years.size()); + canv->Divide(3, years.size()); + + string massCut = Form("(mass>%.2f && mass<%.2f)", + rangeVal[comp->at(0)], + rangeVal[comp->at(0)+1]); + string rangeList = rangeName[comp->at(0)]; + + // cut on mass range + for (uint iCut=1; iCutsize(); ++iCut) { + massCut = massCut + Form("||(mass>%.2f && mass<%.2f)", + rangeVal[comp->at(iCut)], + rangeVal[comp->at(iCut)+1]); + rangeList = rangeList + "," + rangeName[comp->at(iCut)]; + } + + // loop on years + for (unsigned int iy = 0; iy < years.size(); iy++) { + year.clear(); year.assign(Form("%i",years[iy])); + + std::string frametitle = ""; + + std::vector frames; + frametitle = Form("cos#theta_{K} projection of %s dataset, mass range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( ctK ->frame(Title(frametitle.c_str())) )); + frametitle = Form("cos#theta_{L} projection of %s dataset, mass range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( ctL ->frame(Title(frametitle.c_str())) )); + frametitle = Form("#phi projection of %s dataset, mass range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( phi ->frame(Title(frametitle.c_str())) )); + TLegend* leg = new TLegend (0.45,0.7,0.9,0.9); + + auto singleYearPdf = simPdf->getPdf(("data"+year+Form("_subs%d",iSample)).c_str()); + auto singleYearData = combData->reduce(("sample==sample::data"+year+Form("_subs%d",iSample)).c_str()); + + TLatex text_chi2; + text_chi2.SetTextSize(0.04); + + for (unsigned int fr = 0; fr < frames.size(); fr++){ + + singleYearData->plotOn(frames[fr], + MarkerColor(kRed+1), + LineColor(kRed+1), + Binning(40), + Cut(massCut.c_str()), + // Cut(("("+massCut+")&&sample==sample::data"+year+Form("_subs%d",iSample)).c_str()), + Name(("plData"+name+year).c_str())); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + Name(("plPDF"+name+year).c_str()), + NumCPU(8), + RooFit::Precision(1) + ); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(8), + Name(("plPDFbkg"+name+year).c_str()), + Components( ("bkg_pdf_"+year).c_str() ), + NumCPU(8), + RooFit::Precision(1) + ); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(880), + Name(("plPDFsig"+name+year).c_str()), + Components( ("PDF_sig_ang_mass_"+shortString+"_"+year).c_str() ), + NumCPU(8), + RooFit::Precision(1) + ); + + if (fr == 0) { + leg->AddEntry(frames[fr]->findObject(("plData"+name+year).c_str()), + "Data", "lep"); + leg->AddEntry(frames[fr]->findObject(("plPDF"+name+year ).c_str()), + "Total PDF","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFsig"+name+year ).c_str()), + "Signal","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFbkg"+name+year ).c_str()), + "Background","l"); + } + canv->cd(iy*3+fr+1); + gPad->SetLeftMargin(0.16); + frames[fr]->Draw(); + double chiSq = frames[fr]->chiSquare(("plPDF"+name+year).c_str(),("plData"+name+year).c_str()); + text_chi2.DrawLatexNDC(0.22,0.86,Form("#chi^{2}/N_{bins} = %.2f/ 40",chiSq*40)); + if (fr == 0) leg->Draw("same"); + } + } + + canv->SaveAs( ("plotSimFit4d_d/slices/simFitResult_data_" + plotString + "_" + name + "_massSlices.pdf").c_str() ); + + } + + return; + +} + + + +void plot_simfit_data_fullAngularMass_SwaveBin1(int q2Bin, int parity, uint nSample, std::vector years) +{ + if ( parity==-1 ) + for (parity=0; parity<2; ++parity) + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); + else + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); +} + +int main(int argc, char** argv) +{ + // q2-bin format: [0-8] for one bin + // [-1] for each bin recursively + // parity format: [0] even efficiency + // [1] odd efficiency + // [-1] for each parity recursively + + int q2Bin = -1; + int parity = -1; + + if ( argc > 1 ) q2Bin = atoi(argv[1]); + if ( argc > 2 ) parity = atoi(argv[2]); + + uint nSample = 0; + if ( argc > 3 ) nSample = atoi(argv[3]); + + std::vector years; + if ( argc > 4 && atoi(argv[4]) != 0 ) years.push_back(atoi(argv[4])); + else { + cout << "No specific years selected, using default: 2016" << endl; + years.push_back(2016); + } + if ( argc > 5 && atoi(argv[5]) != 0 ) years.push_back(atoi(argv[5])); + if ( argc > 6 && atoi(argv[6]) != 0 ) years.push_back(atoi(argv[6])); + + cout << "q2Bin " << q2Bin << endl; + cout << "parity " << parity << endl; + cout << "nSample " << nSample << endl; + cout << "years[0] " << years[0] << endl; +// cout << years[1] << endl; +// cout << years[2] << endl; + + + if ( q2Bin < -1 || q2Bin >= nBins ) return 1; + if ( parity < -1 || parity > 1 ) return 1; + + if ( q2Bin==-1 ) cout << "Running all the q2 bins" << endl; + if ( parity==-1 ) cout << "Running both the parity datasets" << endl; + + if ( q2Bin==-1 ) + for (q2Bin=0; q2Bin +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "RooDoubleCBFast.h" +#include "RooExponential.h" + +#include "utils.h" +#include "PdfSigRTMass.h" +#include "PdfSigWTMass.h" +#include "PdfSigAngMass.h" +#include "ShapeSigAng.h" + +using namespace RooFit; +using namespace std; + +static const int nBins = 9; + +TCanvas* c [4*nBins]; + +double power = 1.0; + +static const int nRanges = 3; +double rangeVal[nRanges+1] = {-1 * TMath::Pi(), -1 * TMath::Pi()/3, TMath::Pi()/3, TMath::Pi()}; +string rangeName[nRanges] = {"lowPhi","mediumPhi","highPhi"}; + +void plot_simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, uint nSample, std::vector years) +{ + + RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ; + + // to speed up code + RooAbsReal::defaultIntegratorConfig()->getConfigSection("RooIntegrator1D").setRealValue("fixSteps",2) ; + + map< string, vector* > rangeComp; + for (int i=0; i {i}; + + string shortString = Form("b%ip%i",q2Bin,parity); + cout<<"Conf: "<0) stat = Form("_sub-%i",iSample); + string filename_nominal = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/simFitResults/simFitResult_data_fullAngularMass_Swave_201620172018_b%i-XGBv8_unbl4.root",q2Bin); + TFile* fin = new TFile(filename_nominal.c_str(),"READ"); + auto wsp_out = (RooWorkspace*)fin->Get("wsp_out"); + auto combData = (RooDataSet*)wsp_out->data("allcombData"); + auto simPdf = (RooSimultaneous*)wsp_out->pdf("simPdf"); + auto mass = (RooRealVar*)wsp_out->var("mass"); + auto ctK = (RooRealVar*)wsp_out->var("ctK"); + auto ctL = (RooRealVar*)wsp_out->var("ctL"); + auto phi = (RooRealVar*)wsp_out->var("phi"); + + string plotString = shortString + "_" + all_years; + if (nSample>0) plotString = plotString + Form("_s%i",nSample); + + int confIndex = 2*nBins*parity + q2Bin; + string longString = "Fit to reconstructed events"; + longString = longString + Form(parity==1?" (q2-bin %i even)":" (q2-bin %i odd) ",q2Bin); + + + // plot fit projections + c[confIndex] = new TCanvas (("c_"+shortString).c_str(),("Fit to RECO-level MC - "+longString).c_str(),2000,1500); + if (years.size()>1) c[confIndex]->Divide(4, years.size()); + else c[confIndex]->Divide(2,2); + + // now, plot fit projections in different defined ctK ranges + + // assign ranges to the ctK variable + for (int i=0; isetRange(rangeName[i].c_str(),rangeVal[i],rangeVal[i+1]); + + + for (const auto& [name, comp] : rangeComp) { + auto canv = new TCanvas (Form("canv%s",name.c_str()), + Form("canv%s",name.c_str()), + 1500,500*years.size()); + canv->Divide(3, years.size()); + + string phiCut = Form("(phi>%.2f && phi<%.2f)", + rangeVal[comp->at(0)], + rangeVal[comp->at(0)+1]); + string rangeList = rangeName[comp->at(0)]; + + // cut on mass range + for (uint iCut=1; iCutsize(); ++iCut) { + phiCut = phiCut + Form("||(ctK>%.2f && ctK<%.2f)", + rangeVal[comp->at(iCut)], + rangeVal[comp->at(iCut)+1]); + rangeList = rangeList + "," + rangeName[comp->at(iCut)]; + } + + // loop on years + for (unsigned int iy = 0; iy < years.size(); iy++) { + year.clear(); year.assign(Form("%i",years[iy])); + + std::string frametitle = ""; + + std::vector frames; + frametitle = Form("mass projection of %s dataset, #phi range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( mass ->frame(Title(frametitle.c_str())) )); + frametitle = Form("cos#theta_{L} projection of %s dataset, #phi range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( ctL ->frame(Title(frametitle.c_str())) )); + frametitle = Form("cos#theta_{K} projection of %s dataset, #phi range [%.2f,%.2f]",year.c_str(),rangeVal[comp->at(0)], rangeVal[comp->at(0)+1]); + frames.push_back( prepareFrame( ctK ->frame(Title(frametitle.c_str())) )); + TLegend* leg = new TLegend (0.45,0.7,0.9,0.9); + + auto singleYearPdf = simPdf->getPdf(("data"+year+Form("_subs%d",iSample)).c_str()); + auto singleYearData = combData->reduce(("sample==sample::data"+year+Form("_subs%d",iSample)).c_str()); + + TLatex text_chi2; + text_chi2.SetTextSize(0.035); + + for (unsigned int fr = 0; fr < frames.size(); fr++){ + + std::cout << "now working on frame " << fr << std::endl; + singleYearData->plotOn(frames[fr], + MarkerColor(kRed+1), + LineColor(kRed+1), + Binning(40), + Cut(phiCut.c_str()), + // Cut(("("+massCut+")&&sample==sample::data"+year+Form("_subs%d",iSample)).c_str()), + Name(("plData"+name+year).c_str())); + + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + Name(("plPDF"+name+year).c_str()), + NumCPU(8), + RooFit::Precision(1) + ); + + std::cout << "\t full PDF: done" << std::endl; + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(8), + Name(("plPDFbkg"+name+year).c_str()), + Components( ("bkg_pdf_"+year).c_str() ), + NumCPU(8), + RooFit::Precision(1) + ); + + std::cout << "\t bkg PDF: done" << std::endl; + singleYearPdf->plotOn(frames[fr], + // Slice(sample, ("data"+year+Form("_subs%d",iSample)).c_str()), + // ProjWData(RooArgSet(sample), *combData), + ProjectionRange(rangeList.c_str()), + LineWidth(1), + LineColor(880), + Name(("plPDFsig"+name+year).c_str()), + Components( ("PDF_sig_ang_mass_"+shortString+"_"+year).c_str() ), + NumCPU(8), + RooFit::Precision(1) + ); + + std::cout << "\t sig PDF: done" << std::endl; + if (fr == 0) { + leg->AddEntry(frames[fr]->findObject(("plData"+name+year).c_str()), + "Data", "lep"); + leg->AddEntry(frames[fr]->findObject(("plPDF"+name+year ).c_str()), + "Total PDF","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFsig"+name+year ).c_str()), + "Signal","l"); + leg->AddEntry(frames[fr]->findObject(("plPDFbkg"+name+year ).c_str()), + "Background","l"); + } + canv->cd(iy*3+fr+1); + gPad->SetLeftMargin(0.16); + frames[fr]->Draw(); + double chiSq = frames[fr]->chiSquare(("plPDF"+name+year).c_str(),("plData"+name+year).c_str()); + text_chi2.DrawLatexNDC(0.22,0.86,Form("#chi^{2}/N_{bins} = %.2f/ 40",chiSq*40)); + if (fr == 0) leg->Draw("same"); + } + } + + canv->SaveAs( ("plotSimFit4d_d/slices/simFitResult_data_" + plotString + "_" + name + "_phiSlices.pdf").c_str() ); + + } + + return; + +} + + + +void plot_simfit_data_fullAngularMass_SwaveBin1(int q2Bin, int parity, uint nSample, std::vector years) +{ + if ( parity==-1 ) + for (parity=0; parity<2; ++parity) + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); + else + plot_simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, nSample, years); +} + +int main(int argc, char** argv) +{ + // q2-bin format: [0-8] for one bin + // [-1] for each bin recursively + // parity format: [0] even efficiency + // [1] odd efficiency + // [-1] for each parity recursively + + int q2Bin = -1; + int parity = -1; + + if ( argc > 1 ) q2Bin = atoi(argv[1]); + if ( argc > 2 ) parity = atoi(argv[2]); + + uint nSample = 0; + if ( argc > 3 ) nSample = atoi(argv[3]); + + std::vector years; + if ( argc > 4 && atoi(argv[4]) != 0 ) years.push_back(atoi(argv[4])); + else { + cout << "No specific years selected, using default: 2016" << endl; + years.push_back(2016); + } + if ( argc > 5 && atoi(argv[5]) != 0 ) years.push_back(atoi(argv[5])); + if ( argc > 6 && atoi(argv[6]) != 0 ) years.push_back(atoi(argv[6])); + + cout << "q2Bin " << q2Bin << endl; + cout << "parity " << parity << endl; + cout << "nSample " << nSample << endl; + cout << "years[0] " << years[0] << endl; +// cout << years[1] << endl; +// cout << years[2] << endl; + + + if ( q2Bin < -1 || q2Bin >= nBins ) return 1; + if ( parity < -1 || parity > 1 ) return 1; + + if ( q2Bin==-1 ) cout << "Running all the q2 bins" << endl; + if ( parity==-1 ) cout << "Running both the parity datasets" << endl; + + if ( q2Bin==-1 ) + for (q2Bin=0; q2Bin0) fileName = fileDir + "simFitResults4d/simFitResult_recoMC_fullAngularMass" + all_years + Form("_dataStat-%i_b%i.root", nSample, q2Bin); TFile* fin = new TFile(fileName.c_str(),"READ"); auto wsp_out = (RooWorkspace*)fin->Get("wsp_out"); diff --git a/plot_utils.h b/plot_utils.h new file mode 100644 index 0000000..bb7a4e1 --- /dev/null +++ b/plot_utils.h @@ -0,0 +1,169 @@ +// common definitions and functions + +static const int nBins = 8; // was 8 +float binBorders [nBins+1] = { 1.1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16}; +// static const int nBins = 9; +// float binBorders [nBins+1] = { 1.1, 2, 4.3, 6, 8.68, 10.09, 12.86, 14.18, 16, 19}; + +static const int nPars = 8; +string ParName [nPars] = { "Fl", "P1", "P2", "P3", "P4p", "P5p", "P6p", "P8p" }; +string ParTitle [nPars] = { "F_{L}", "P_{1}", "P_{2}", "P_{3}", "P'_{4}", "P'_{5}", "P'_{6}", "P'_{8}" }; +double ParMin [nPars] = {0,-1,-0.5,-0.5,-sqrt(2),-sqrt(2),-sqrt(2),-sqrt(2)}; +double ParMax [nPars] = {1,1,0.5,0.5,sqrt(2),sqrt(2),sqrt(2),sqrt(2)}; + + +void setGraphProperties(TGraphAsymmErrors* Gr, int icolor, int thiscolor, string label, int ParIndx, string type, int marker =20){ + + int addc = 0; + if (type.find("WT") != std::string::npos) addc = 1; + if (type.find("CT") != std::string::npos) addc = 2; + + Gr -> SetLineColor(thiscolor + addc); + Gr -> SetMarkerColor(thiscolor + addc); +// Gr -> SetLineColor(colorlist[icolor] + addc); +// Gr -> SetMarkerColor(colorlist[icolor] + addc); + Gr -> SetMarkerStyle(marker + icolor + addc); + Gr -> SetName(Form("Gr%s%i_%s" , type.c_str(), ParIndx, label.c_str()) ); + Gr -> SetLineWidth(2); +} + +void setAxisProperties(TH1F* auxE2, float diffMax){ + auxE2->SetStats(kFALSE); + auxE2->SetLineColor(1); + auxE2->GetXaxis()->SetTitle("q^{2} (GeV^{2})"); + auxE2->GetXaxis()->SetTitleSize(0.12); + auxE2->GetXaxis()->SetTitleOffset(0.95); + auxE2->GetXaxis()->SetLabelSize( 0.10); + auxE2->GetXaxis()->SetTickLength(0.1); + auxE2->GetYaxis()->SetTitle("(RECO - GEN)"); + auxE2->GetYaxis()->SetTitleSize(0.10); + auxE2->GetYaxis()->SetTitleOffset(0.45); + auxE2->GetYaxis()->SetRangeUser(-1*diffMax,diffMax); + auxE2->GetYaxis()->SetLabelSize(0.09); + auxE2->GetYaxis()->SetNdivisions(505); +} + +TGraphErrors* getGreyResonanceGraph(int ParIndx){ + + double ResX [2] = {0.5*(binBorders[5]+binBorders[4]),0.5*(binBorders[7]+binBorders[6])}; + double ResXe[2] = {0.5*(binBorders[5]-binBorders[4]+0.01),0.5*(binBorders[7]-binBorders[6])}; + double ResY [2] = {0.5*(ParMax[ParIndx]+ParMin[ParIndx])-0.002*(ParMax[ParIndx]-ParMin[ParIndx]), + 0.5*(ParMax[ParIndx]+ParMin[ParIndx])-0.002*(ParMax[ParIndx]-ParMin[ParIndx])}; + double ResYe[2] = {0.498*(ParMax[ParIndx]-ParMin[ParIndx]),0.498*(ParMax[ParIndx]-ParMin[ParIndx])}; + TGraphErrors *resCover = new TGraphErrors(2,ResX,ResY,ResXe,ResYe); + resCover->SetName(Form("resCover%i",ParIndx)); + resCover->SetFillColor(18); + resCover->SetFillStyle(1001); + return resCover; +} + +TGraphErrors* getGreyResonancePull(int ParIndx, float diffMax){ + + double ResX [2] = {0.5*(binBorders[5]+binBorders[4]),0.5*(binBorders[7]+binBorders[6])}; + double ResXe[2] = {0.5*(binBorders[5]-binBorders[4]+0.01),0.5*(binBorders[7]-binBorders[6])}; + double ResD [2] = {-0.01*diffMax,-0.01*diffMax}; + double ResDe[2] = { 0.99*diffMax, 0.99*diffMax}; + TGraphErrors *resDiffCover = new TGraphErrors(2,ResX,ResD,ResXe,ResDe); + resDiffCover->SetName(Form("resDiffCover%i",ParIndx)); + resDiffCover->SetFillColor(18); + resDiffCover->SetFillStyle(1001); + return resDiffCover; +} + +void drawBinLines (float diffMax){ + std::vector lines; + for (int i=0; iSetLineStyle(3); + lines[i]->SetLineColor(kGray); + lines[i]->Draw(); + } + return; +} +// fill arrays with fit results, reading from RooFitResult +// compute difference wrt reference values (genRes here) +void fillVectors(RooFitResult* fitResult, int ParIndx, double Res[][nBins], double ErrH[][nBins], double ErrL[][nBins], + double Diff[][nBins], double DErrH[][nBins], double DErrL[][nBins], + double genRes[nBins], double genErrH[nBins], double genErrL[nBins], + int fileN, int ibin){ + RooRealVar* Par = (RooRealVar*)fitResult->floatParsFinal().find(ParName[ParIndx].c_str()); + Res[fileN][ibin] = Par->getValV(); + ErrH[fileN][ibin] = Par->getErrorHi(); + ErrL[fileN][ibin] = fabs(Par->getErrorLo()); // return signed value in GEN fit results + Diff[fileN][ibin] = Res[fileN][ibin]-genRes[ibin]; + DErrH[fileN][ibin] = sqrt(ErrH[fileN][ibin]*ErrH[fileN][ibin]+genErrL[ibin]*genErrL[ibin]); + DErrL[fileN][ibin] = sqrt(ErrL[fileN][ibin]*ErrL[fileN][ibin]+genErrH[ibin]*genErrH[ibin]); +} + +// fill arrays with fit results, reading from RooFitResult +void fillVectors(RooFitResult* fitResult, int ParIndx, double Res[nBins], double ErrH[nBins], double ErrL[nBins], int ibin){ + RooRealVar* Par = (RooRealVar*)fitResult->floatParsFinal().find(ParName[ParIndx].c_str()); + Res[ibin] = Par->getValV(); + ErrH[ibin] = Par->getErrorHi(); + ErrL[ibin] = fabs(Par->getErrorLo()); +// ErrL[ibin] = -1*Par->getErrorLo(); +// std::cout << Par->GetName() << " " << ibin << " " << Res[ibin] << " + " << ErrH[ibin] << " - " << ErrL[ibin]<< std::endl; +} + + +// fill arrays with fit results, reading from TTree containing fitted values +// compute difference wrt reference values (genRes here) +void fillVectors(TTree* fr, int ParIndx, + double Res[][nBins], double ErrH[][nBins], double ErrL[][nBins], + double Diff[][nBins], double DErrH[][nBins], double DErrL[][nBins], + double genRes[nBins], double genErrH[nBins], double genErrL[nBins], + int fileN, int ibin, int nentry = 0, int debug = false){ + + double val{100}; + double ci[2]; // confidence interval extremes + fr->SetBranchAddress(Form("%s_best",ParName[ParIndx].c_str()),&val); + fr->SetBranchAddress(Form("%s_high",ParName[ParIndx].c_str()),&ci[0]); + fr->SetBranchAddress(Form("%s_low" ,ParName[ParIndx].c_str()),&ci[1]); + + fr->GetEntry(0); + Res[fileN][ibin] = val; + ErrH[fileN][ibin] = ci[0]-val; + ErrL[fileN][ibin] = val - ci[1]; + if (nentry==-1){ + ErrH[fileN][ibin] = ci[0]; + ErrL[fileN][ibin] = fabs(ci[1]); + } + if (debug){ + std::cout << ParName[ParIndx].c_str() << " bin: " << ibin << " : " << val << " + " << ErrH[fileN][ibin] << " - " << ErrL[fileN][ibin]<< std::endl; + std::cout << "\t diff to GEN " << Res[fileN][ibin]-genRes[ibin]<< std::endl; + } + Diff[fileN][ibin] = Res[fileN][ibin]-genRes[ibin]; + DErrH[fileN][ibin] = sqrt(ErrH[fileN][ibin]*ErrH[fileN][ibin]+genErrL[ibin]*genErrL[ibin]); + DErrL[fileN][ibin] = sqrt(ErrL[fileN][ibin]*ErrL[fileN][ibin]+genErrH[ibin]*genErrH[ibin]); +} + + +// fill arrays with fit results (average and RMS), reading from TChain containing N fitted values +// compute difference wrt reference values (genRes here) +void fillVectors(TChain* fr, int ParIndx, + double Res[][nBins], double ErrH[][nBins], double ErrL[][nBins], + double Diff[][nBins], double DErrH[][nBins], double DErrL[][nBins], + double genRes[nBins], double genErrH[nBins], double genErrL[nBins], + int fileN, int ibin, int nentry = 0){ + + double val{100}; + double ci[2]; // confidence interval extremes + fr->SetBranchAddress(Form("%s_best",ParName[ParIndx].c_str()),&val); + fr->SetBranchAddress(Form("%s_high",ParName[ParIndx].c_str()),&ci[0]); + fr->SetBranchAddress(Form("%s_low" ,ParName[ParIndx].c_str()),&ci[1]); + + fr->Draw(Form("%s_best",ParName[ParIndx].c_str())); + TH1F *htemp = (TH1F*)gPad->GetPrimitive("htemp"); + Res[fileN][ibin] = htemp->GetMean(); + ErrH[fileN][ibin] = htemp->GetRMS();///fr->GetEntries(); + ErrL[fileN][ibin] = htemp->GetRMS();///fr->GetEntries(); + + Diff[fileN][ibin] = Res[fileN][ibin]-genRes[ibin]; + DErrH[fileN][ibin] = sqrt(ErrH[fileN][ibin]*ErrH[fileN][ibin]+genErrL[ibin]*genErrL[ibin]); + DErrL[fileN][ibin] = sqrt(ErrL[fileN][ibin]*ErrL[fileN][ibin]+genErrH[ibin]*genErrH[ibin]); + +} + + + + diff --git a/run_simfit_data_fullAngularMass_Swave.sh b/run_simfit_data_fullAngularMass_Swave.sh index abd78c7..344e86c 100755 --- a/run_simfit_data_fullAngularMass_Swave.sh +++ b/run_simfit_data_fullAngularMass_Swave.sh @@ -1,22 +1,19 @@ #!/bin/bash - -par=1 - multi=0 -nsam=${1} -q2stat=${4} - XGBv=8 -localFile=0 -fitopt=0 plot=1 -save=${5} -bin=${2} -# ibin=${2} +localFile=0 +fitopt=0 +unbl=4 +nsam=${1} +bin=${2} yearConf=${3} +q2stat=${4} +save=${5} +par=${6} export SAMPLEDIR=/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/data-datasets/ export SBDIR=/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/sidebands/ @@ -102,26 +99,26 @@ mkdir simFitResults4d mkdir plotSimFit4d_d case "$yearConf" in -# void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSample, uint nSample, uint q2stat, int fitOption, int XGBv, bool localFiles, bool plot, int save, std::vector years) +# void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSample, uint nSample, uint q2stat, int fitOption, int XGBv, int unblind, bool localFiles, bool plot, int save, std::vector years) 0) - echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2016 2017 2018 - ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2016 2017 2018 + echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 2017 2018 + ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 2017 2018 ;; 1) - echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2016 - ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2016 + echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 + ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 ;; 2) - echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2017 - ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2017 + echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2017 + ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2017 ;; 3) - echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2018 - ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${localFile} ${plot} ${save} 2018 + echo ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2018 + ./simfit_data_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2018 ;; esac diff --git a/run_simfit_recoMC_fullAngular.sh b/run_simfit_recoMC_fullAngular.sh index 15a09e7..89921d5 100755 --- a/run_simfit_recoMC_fullAngular.sh +++ b/run_simfit_recoMC_fullAngular.sh @@ -1,16 +1,14 @@ #!/bin/bash - -par=1 - multi=0 -nsam=${1} - xgb=8 plot=0 save=1 +nsam=${1} ibin=${2} +par=${3} +yearConf=${4} ## if localFile == 1, please check the name of the MC dataset being copied localFile=0 @@ -96,8 +94,28 @@ cp $HOME/simfit_recoMC_fullAngular . mkdir -p simFitResults/xgbv8 mkdir -p plotSimFit_d/xgbv8 -echo ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 -./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 +case "$yearConf" in + 0) + echo ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 + ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 + ;; + + 1) + echo ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 + ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 + ;; + + 2) + echo ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2017 + ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2017 + ;; + + 3) + echo ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2018 + ./simfit_recoMC_fullAngular ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2018 + ;; + +esac if [ ! -d $HOME/simFitResults/xgbv8 ]; then mkdir -p $HOME/simFitResults/xgbv8 @@ -107,6 +125,7 @@ if [ ! -d $HOME/plotSimFit_d/xgbv8 ]; then fi cp plotSimFit_d/xgbv8/* $HOME/plotSimFit_d/xgbv8/ cp simFitResults/xgbv8/* $HOME/simFitResults/xgbv8/ + # for file in simFitResults/* ; do cp $file $HOME/${file//.root/_${multi}s${nsam}.root}; done # rm -rf plotSimFit_d diff --git a/run_simfit_recoMC_fullAngularMass.sh b/run_simfit_recoMC_fullAngularMass.sh index 8c5f8f9..9aebbd0 100755 --- a/run_simfit_recoMC_fullAngularMass.sh +++ b/run_simfit_recoMC_fullAngularMass.sh @@ -1,18 +1,17 @@ #!/bin/bash - -par=1 - multi=0 -nsam=${1} - xgb=8 plot=1 save=1 +nsam=${1} ibin=${2} +par=${3} +yearConf=${4} localFile=0 + export SAMPLEDIR=/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/MC-datasets/ if [ "${xgb}" == 0 ]; then @@ -93,14 +92,28 @@ cp $HOME/*.pcm . mkdir -p simFitResults4d/xgbv8 mkdir -p plotSimFit4d_d/xgbv8 -if [ "$bin" -lt 8 ]; then - echo ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 - ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 -else - echo ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 - ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 -fi +case "$yearConf" in + 0) + echo ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 + ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 2017 2018 + ;; + + 1) + echo ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 + ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2016 + ;; + + 2) + echo ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2017 + ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2017 + ;; + + 3) + echo ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2018 + ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${xgb} ${localFile} ${plot} ${save} 2018 + ;; +esac if [ ! -d $HOME/simFitResults4d/xgbv8/ ]; then mkdir -p $HOME/simFitResults4d/xgbv8/ diff --git a/run_simfit_recoMC_fullAngularMass_toybkg.sh b/run_simfit_recoMC_fullAngularMass_toybkg.sh index 45a9d3c..13907af 100755 --- a/run_simfit_recoMC_fullAngularMass_toybkg.sh +++ b/run_simfit_recoMC_fullAngularMass_toybkg.sh @@ -1,16 +1,13 @@ #!/bin/bash - -par=1 - multi=0 -nsam=${1} - xgb=8 plot=1 save=1 +nsam=${1} ibin=${2} +par=${3} ## if localFile == 1, please check the name of the MC dataset being copied localFile=0 diff --git a/run_simfit_toy_fullAngularMass_Swave.sh b/run_simfit_toy_fullAngularMass_Swave.sh new file mode 100644 index 0000000..2a42d99 --- /dev/null +++ b/run_simfit_toy_fullAngularMass_Swave.sh @@ -0,0 +1,137 @@ +#!/bin/bash + +par=0 + +multi=0 +nsam=${1} +q2stat=0 + +XGBv=8 +localFile=0 +fitopt=0 +unbl=4 + +plot=0 +save=1 + +bin=${2} +# ibin=${2} + +yearConf=0 + +export SAMPLEDIR=/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/data-datasets/ +export SBDIR=/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/sidebands/ + +if [ "${XGBv}" == 0 ]; then + export EFFDIR=/eos/user/a/aboletti/BdToKstarMuMu/eff-KDE-theta-v7/files +else + export EFFDIR=/eos/user/a/aboletti/BdToKstarMuMu/eff-KDE-theta-v7-XGBv8/files +fi + +if [ "${USER}" == "aboletti" ]; then + export HOME=/afs/cern.ch/work/a/aboletti/private/Kstmumu-Run2/UML-fit-March2024 + export CMSSWDIR=/afs/cern.ch/work/a/aboletti/private/Kstmumu-Run2/CMSSW_10_4_0/src +elif [ "${USER}" == "fiorendi" ]; then + export HOME=/afs/cern.ch/work/f/fiorendi/private/effKDE/UML-fit + export CMSSWDIR=/afs/cern.ch/work/f/fiorendi/private/effKDE/CMSSW_10_4_0/src +else + echo no user found + exit 1 +fi + +export WORKDIR=$PWD +cd $CMSSWDIR +source /cvmfs/cms.cern.ch/cmsset_default.sh +eval `scram runtime -sh` + +echo setting HOME to $HOME +echo setting CMSSWDIR to $CMSSWDIR + +cd $WORKDIR + +# nbin=0 +# while read -a line; do +# abin[$nbin]=${line[0]} +# nbin=$((nbin+1)) +# done < $HOME/../confSF/KDE_SF.list +# bin=${abin[$ibin]} + +echo 'now submitting for bin ' ${bin} + +if [ "${par}" == 0 ]; then + parstr="ev" +else + parstr="od" +fi + +if [ "${localFile}" -gt 0 ]; then + for iy in {2016..2018} + do + echo 'will copy data samples from from ' ${SAMPLEDIR} + echo 'will copy efficiencies from ' ${EFFDIR} + echo 'will copy sidebands from ' ${SBDIR} + [ "$yearConf" -gt 0 ] && [ "$((${yearConf}+2015))" != "$iy" ] && continue + dataname="${SAMPLEDIR}/recoDATADataset_b${bin}_${iy}.root" + effname="${EFFDIR}/KDEeff_b${bin}_${parstr}_${iy}.root" + sbname="${SBDIR}/savesb_${iy}_b${bin}.root" + if [ ! -r "${dataname}" ]; then + echo "${dataname}" not found + exit 1 + fi + if [ ! -r "${effname}" ]; then + echo "${effname}" not found + exit 1 + fi + if [ ! -r "${sbname}" ]; then + echo "${sbname}" not found + exit 1 + fi + cp "${dataname}" . + cp "${effname}" . + cp "${sbname}" . + done +fi + +if [ ! -r $HOME/simfit_toy_fullAngularMass_Swave ]; then + echo $HOME/simfit_toy_fullAngularMass_Swave not found + exit 1 +fi +cp $HOME/simfit_toy_fullAngularMass_Swave . +cp $HOME/*.pcm . + +mkdir simFitResults4d +mkdir plotSimFit4d_d + +case "$yearConf" in +# void simfit_toy_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSample, uint nSample, uint q2stat, int fitOption, int XGBv, int unblind, bool localFiles, bool plot, int save, std::vector years) + + 0) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 2017 2018 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 2017 2018 + ;; + + 1) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 + ;; + + 2) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2017 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2017 + ;; + + 3) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2018 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2018 + ;; + +esac + +if [ ! -d $HOME/simFitResults4d ]; then + mkdir $HOME/simFitResults4d +fi +if [ ! -d $HOME/plotSimFit4d_d ]; then + mkdir $HOME/plotSimFit4d_d +fi +cp plotSimFit4d_d/* $HOME/plotSimFit4d_d/ +cp simFitResults4d/* $HOME/simFitResults4d/ diff --git a/simfit_data_fullAngularMass_Swave.cc b/simfit_data_fullAngularMass_Swave.cc index 61cacec..f2e2107 100644 --- a/simfit_data_fullAngularMass_Swave.cc +++ b/simfit_data_fullAngularMass_Swave.cc @@ -56,6 +56,8 @@ TCanvas* cPen; TCanvas* c [4*nBins]; double power = 1.0; +bool evalMistagSys = false; + void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSample, uint nSample, uint q2stat, int fitOption, int XGBv, int unblind, bool localFiles, bool plot, int save, std::vector years) { @@ -240,6 +242,8 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl if ( !effCHist[iy] || effCHist[iy]->IsZombie() || !effWHist[iy] || effWHist[iy]->IsZombie() ) { cout<<"Efficiency histogram "<< effCString <<" or " << effWString << " not found in file: "<< filename < R = 1 +/- " << frac_sigma << endl; + cout << "mFrac = " << fraction << " +/- " << fM_sigmas[years[iy]][q2Bin] << " ---> R = " << mFrac->getVal() << " +/- " << frac_sigma << endl; c_vars.add(*mFrac); // Angular Component @@ -521,6 +532,7 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl // in Jpsi bin, if fitOption > 0 -> bernstein polynbomial RooRealVar* slope = new RooRealVar(Form("slope^{%i}",years[iy]), Form("slope^{%i}",years[iy]) , -5., -10., 0.); RooAbsPdf* bkg_mass_pdf = 0; + RooExponential* bkg_exp_pdf = 0; double pol_bmax =1.; RooRealVar* b0_bkg_mass = new RooRealVar(Form("b0_bkg_mass-%i",years[iy]) , Form("b0_bkg_mass-%i",years[iy]) , pol_bmax ); RooRealVar* b1_bkg_mass = new RooRealVar(Form("b1_bkg_mass-%i",years[iy]) , Form("b1_bkg_mass-%i",years[iy]) , 0.1, 0., pol_bmax); @@ -532,7 +544,13 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl b3_bkg_mass->setConstant(true); bkg_mass_pdf = new RooBernstein(Form("bkg_mass1_%i",years[iy]),Form("bkg_mass1_%i",years[iy]), *mass, RooArgList(*b0_bkg_mass, *b1_bkg_mass, *b2_bkg_mass, *b3_bkg_mass,* b4_bkg_mass)); } else { - bkg_mass_pdf = new RooExponential(Form("bkg_mass1_%i",years[iy]), Form("bkg_mass1_%i",years[iy]) , *mass, *slope ); + bkg_exp_pdf = new RooExponential(Form("bkg_mass1_%i",years[iy]), Form("bkg_mass1_%i",years[iy]) , *mass, *slope ); + if (q2Bin==3||q2Bin==5) { + bkg_mass_pdf = createPdfCuts(q2Bin,years[iy], mass, slope); + }else{ + bkg_mass_pdf = bkg_exp_pdf; + } + } RooProdPdf* bkg_pdf = new RooProdPdf(Form(bkgpdfname.c_str(),years[iy]), @@ -627,7 +645,9 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl } TFile* fout = 0; - if (save>0) fout = new TFile(("simFitResults4d/simFitResult_data_fullAngularMass_Swave_" + all_years + stat + Form("_b%i%s_unbl%i.root", q2Bin, XGBstr.c_str(), unblind)).c_str(),"RECREATE"); + std::string addInfo = evalMistagSys ? "_mistagSyst_exp" : ""; + + if (save>0) fout = new TFile(("simFitResults4d/simFitResult_data_fullAngularMass_Swave_" + all_years + stat + Form("_b%i%s_unbl%i%s.root", q2Bin, XGBstr.c_str(), unblind, addInfo.c_str())).c_str(),"RECREATE"); RooWorkspace* wsp_out = 0; // save initial par values @@ -647,6 +667,7 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl bool boundCheck, convCheck, usedPenalty; RooArgList pars (*Fl,*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p); + RooArgList sPars (*Fs,*As,*A4s,*A5s,*A7s,*A8s); // TTree with the MINOS output vector vResult (pars.getSize()); @@ -711,7 +732,7 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl *params = *savedParams ; // run the fit - fitter = new Fitter (Form("fitter%i",is),Form("fitter%i",is),pars,combData,simPdf,simPdf_penalty,boundary,bound_dist,penTerm,&c_vars); + fitter = new Fitter (Form("fitter%i",is),Form("fitter%i",is),pars,sPars,combData,simPdf,simPdf_penalty,boundary,bound_dist,penTerm,&c_vars); vFitter.push_back(fitter); fitter->setUnbl(unblind); @@ -725,6 +746,11 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl if (q2Bin==4) fitter->runSimpleFit = true; } +// if (q2Bin==3 || q2Bin==5){ +// fitter->DisableConstOpt(); +// std::cout<fit(); subTime.Stop(); @@ -762,7 +788,7 @@ void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSampl TStopwatch improvTime; improvTime.Start(true); - if (runPostFitSteps) fitter->improveAng(); + if (runPostFitSteps) fitter->improveAng(1,300000); improvTime.Stop(); imprTime = improvTime.CpuTime(); cout<<"Improv time: "<2) { string plotString = shortString + "_" + all_years + stat + XGBstr + Form("_unbl%i",unblind); - string plotname = "plotSimFit4d_d/simFitResult_data_fullAngularMass_Swave_" + plotString + ".pdf"; + addInfo = evalMistagSys ? "_mistagSyst_exp" : ""; + string plotname = "plotSimFit4d_d/simFitResult_data_fullAngularMass_Swave_" + plotString + addInfo + ".pdf"; fitter->plotSimFitProjections(plotname.c_str(),{samplename,sigpdfname,bkgpdfname},years,true); } diff --git a/simfit_data_fullAngularMass_Swave.sh b/simfit_data_fullAngularMass_Swave.sh index eafc17e..9620306 100755 --- a/simfit_data_fullAngularMass_Swave.sh +++ b/simfit_data_fullAngularMass_Swave.sh @@ -1,6 +1,6 @@ #!/bin/bash -par=1 +par=0 multi=0 nsam=0 diff --git a/simfit_recoMC_fullAngularMass.cc b/simfit_recoMC_fullAngularMass.cc index 857985e..72427ab 100644 --- a/simfit_recoMC_fullAngularMass.cc +++ b/simfit_recoMC_fullAngularMass.cc @@ -451,7 +451,7 @@ void simfit_recoMC_fullAngularMassBin(int q2Bin, int parity, bool multiSample, u if (nSample>0) stat = stat + Form("-%i",firstSample); if (multiSample) stat = stat + Form("-%i",lastSample); TFile* fout = 0; - if (save>0) fout = new TFile(("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass" + all_years + stat + Form("_b%i", q2Bin) + XGBstr + ".root").c_str(),"RECREATE"); + if (save>0) fout = new TFile(("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass" + all_years + stat + "_" + shortString + XGBstr + ".root").c_str(),"RECREATE"); RooWorkspace* wsp_out = 0; // save initial par values diff --git a/simfit_recoMC_fullAngularMass.sh b/simfit_recoMC_fullAngularMass.sh index 1f967ff..13e72e0 100755 --- a/simfit_recoMC_fullAngularMass.sh +++ b/simfit_recoMC_fullAngularMass.sh @@ -12,6 +12,8 @@ nsam=0 plot=1 save=2 +XGBv=8 + # Create directories for fit logs, results and plots if [ ! -d logs_simFit4d ]; then mkdir logs_simFit4d; fi if [ ! -d simFitResults4d ]; then mkdir simFitResults4d; fi @@ -24,15 +26,15 @@ if make simfit_recoMC_fullAngularMass; then while read -a line; do bin=${line[0]} - # for year in {2016..2018}; do - - # ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} 0 ${plot} ${save} ${year} \ - # &>logs_simFit4d/simfit_recoMC_fullAngularMass_${bin}_${par}_${multi}_${nsam}_${year}.out & - - # done + for year in {2016..2018}; do + + ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${XGBv} 0 ${plot} ${save} ${year} \ + &>logs_simFit4d/simfit_recoMC_fullAngularMass_${bin}_${par}_${multi}_${nsam}_${XGBv}_${year}.out & + + done - ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} 0 ${plot} ${save} 2016 2017 2018 \ - &>logs_simFit4d/simfit_recoMC_fullAngularMass_${bin}_${par}_${multi}_${nsam}_2016_2017_2018.out & + ./simfit_recoMC_fullAngularMass ${bin} ${par} ${multi} ${nsam} ${XGBv} 0 ${plot} ${save} 2016 2017 2018 \ + &>logs_simFit4d/simfit_recoMC_fullAngularMass_${bin}_${par}_${multi}_${nsam}_${XGBv}_2016_2017_2018.out & done < ../confSF/KDE_SF.list diff --git a/simfit_recoMC_fullAngularMass_toybkg.cc b/simfit_recoMC_fullAngularMass_toybkg.cc index 44cb376..82126e0 100644 --- a/simfit_recoMC_fullAngularMass_toybkg.cc +++ b/simfit_recoMC_fullAngularMass_toybkg.cc @@ -261,10 +261,16 @@ void simfit_recoMC_fullAngularMassBin(int q2Bin, int parity, bool multiSample, u // read mass pdf for background RooRealVar* slope = new RooRealVar (Form("slope^{%i}",years[iy]), Form("slope^{%i}",years[iy]) , wsp_sb[iy]->var("slope")->getVal(), -10., 0.); - RooExponential* bkg_exp = new RooExponential(Form("bkg_exp_%i",years[iy]), Form("bkg_exp_%i",years[iy]) , *slope, *mass ); cout << Form("exponential slope for %f", slope->getVal()) << endl; - - + RooAbsPdf* bkg_mass_pdf = 0; + + RooExponential* bkg_exp_pdf = new RooExponential(Form("bkg_mass1_%i",years[iy]), Form("bkg_mass1_%i",years[iy]) , *mass, *slope ); + if (q2Bin==3 || q2Bin==5) { + bkg_mass_pdf = createPdfCuts(q2Bin,years[iy], mass, slope); + }else{ + bkg_mass_pdf = bkg_exp_pdf; + } + // retrieve sideband range from input file float max_lsb = wsp_sb[iy]->var(Form("max_sbl_bin%i_%i", q2Bin, years[iy]))->getVal(); float min_rsb = wsp_sb[iy]->var(Form("min_sbr_bin%i_%i", q2Bin, years[iy]))->getVal(); @@ -273,8 +279,11 @@ void simfit_recoMC_fullAngularMassBin(int q2Bin, int parity, bool multiSample, u cout << Form("sideband mass range: [5., %.2f] U [%.2f, 5.6]", max_lsb, min_rsb) << endl; // create 4D pdf for background and import to workspace - RooProdPdf* bkg_pdf = new RooProdPdf(Form(bkgpdfname.c_str(),years[iy]), Form(bkgpdfname.c_str(),years[iy]), - RooArgList(*bkg_ang_pdf,*bkg_exp)); + RooProdPdf* bkg_pdf = new RooProdPdf(Form(bkgpdfname.c_str(),years[iy]), + Form(bkgpdfname.c_str(),years[iy]), + RooArgList(*bkg_ang_pdf,*bkg_mass_pdf)); + + wksp->import(*bkg_ang_pdf); wksp->import(*bkg_pdf, RecycleConflictNodes()); @@ -584,7 +593,7 @@ void simfit_recoMC_fullAngularMassBin(int q2Bin, int parity, bool multiSample, u if (nSample>0) stat = stat + Form("-%i",firstSample); if (multiSample) stat = stat + Form("-%i",lastSample); TFile* fout = 0; - if (save>0) fout = new TFile(("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass_toybkg" + all_years + stat + Form("_p%i_b%i.root", parity, q2Bin)).c_str(),"RECREATE"); + if (save>0) fout = new TFile(("simFitResults4d/xgbv8/simFitResult_recoMC_fullAngularMass_toybkg" + all_years + stat + Form("_p%i_b%i", parity, q2Bin) + XGBstr + ".root" ).c_str(),"RECREATE"); RooWorkspace* wsp_out = 0; wksp->import(*simPdf,RecycleConflictNodes()); diff --git a/simfit_toy_fullAngularMass_Swave.cc b/simfit_toy_fullAngularMass_Swave.cc new file mode 100644 index 0000000..86404d1 --- /dev/null +++ b/simfit_toy_fullAngularMass_Swave.cc @@ -0,0 +1,990 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "RooDoubleCBFast.h" +#include "RooExponential.h" +#include "RooPolynomial.h" +#include "RooGenericPdf.h" +#include "RooBernstein.h" + +#include "utils.h" +#include "PdfSigRTMass.h" +#include "PdfSigWTMass.h" +#include "PdfSigAngMass.h" +#include "ShapeSigAng.h" + +#include "BoundCheck.h" +#include "BoundDist.h" +#include "Penalty.h" +#include "Fitter.h" +#include "RooBernsteinSideband.h" + +using namespace RooFit; +using namespace std; + +static const int nBins = 9; + +TCanvas* cnll; +TCanvas* cZoom; +TCanvas* cPen; +TCanvas* c [4*nBins]; + +double power = 1.0; +bool evalMistagSys = false; + + +void simfit_data_fullAngularMass_SwaveBin(int q2Bin, int parity, bool multiSample, uint nSample, uint q2stat, int fitOption, int XGBv, int unblind, bool localFiles, bool plot, int save, std::vector years) +{ + + RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ; + + string shortString = Form("b%ip%i",q2Bin,parity); + cout<<"Conf: "<0) stat = Form("_s%i",nSample); + + string sigpdfname = "PDF_sig_ang_mass_"+shortString+"_%i"; + string bkgpdfname = "bkg_pdf_%i"; + + string XGBstr = ""; + if (XGBv>9) XGBstr = Form("-TMVAv%i",XGBv-10); + else if (XGBv>0) XGBstr = Form("-XGBv%i",XGBv); + + std::vector fin_eff; + std::vector wsp, wsp_mcmass, wsp_sb, wsp_Z4430; + std::vector> data; + std::vector effC, effW; + std::vector effCHist, effWHist; + std::vector intCHist, intWHist; + std::vector< std::vector > intCVec(years.size(), std::vector(0)); + std::vector< std::vector > intWVec(years.size(), std::vector(0)); + std::vector PDF_sig_ang_mass_bkg(0); + std::vector PDF_sig_ang_mass_bkg_penalty(0); + std::vector c_deltaPeaks, c_fm; + RooArgSet c_vars_rt, c_pdfs_rt; + RooArgSet c_vars_wt, c_pdfs_wt; + RooArgSet c_vars; + + //// from https://root-forum.cern.ch/t/combining-roodatasets-using-std-map-in-pyroot/16471/20 + gInterpreter->GenerateDictionary("std::pair", "map;string;RooDataSet.h"); + gInterpreter->GenerateDictionary("std::map", "map;string;RooDataSet.h"); + gInterpreter->GenerateDictionary("std::pair::iterator, bool>", "map;string;RooDataSet.h"); + std::map map; + + RooRealVar* ctK = new RooRealVar("ctK", "cos(#theta_{K})", -1 , 1 ); + RooRealVar* ctL = new RooRealVar("ctL", "cos(#theta_{l})", -1 , 1 ); + RooRealVar* phi = new RooRealVar("phi", "#phi", -3.14159, 3.14159 ); + RooArgList vars (* ctK,* ctL,* phi); + RooRealVar* mass = new RooRealVar("mass","mass", 5.,5.6,"GeV"); + RooRealVar* rand = new RooRealVar("rand", "rand", 0,1); + RooArgSet reco_vars (*ctK, *ctL, *phi, *mass, *rand); + RooArgSet observables (*ctK, *ctL, *phi, *mass); + + // define angular parameters with ranges from positiveness requirements on the decay rate + RooRealVar* Fl = new RooRealVar("Fl","F_{L}",0.5,0,1); + RooRealVar* P1 = new RooRealVar("P1","P_{1}",0,-1,1); + RooRealVar* P2 = new RooRealVar("P2","P_{2}",0,-0.5,0.5); + RooRealVar* P3 = new RooRealVar("P3","P_{3}",0,-0.5,0.5); + RooRealVar* P4p = new RooRealVar("P4p","P'_{4}",0,-1*sqrt(2),sqrt(2)); + RooRealVar* P5p = new RooRealVar("P5p","P'_{5}",0,-1*sqrt(2),sqrt(2)); + RooRealVar* P6p = new RooRealVar("P6p","P'_{6}",0,-1*sqrt(2),sqrt(2)); + RooRealVar* P8p = new RooRealVar("P8p","P'_{8}",0,-1*sqrt(2),sqrt(2)); + + RooRealVar* Fs = new RooRealVar("Fs","F_{S}",0.05,0,1); + RooRealVar* As = new RooRealVar("As","A^{S}",0,-1,1); + RooRealVar* A4s = new RooRealVar("A4s","A_{4}^{S}",0,-1,1); + RooRealVar* A5s = new RooRealVar("A5s","A_{5}^{S}",0,-1,1); + RooRealVar* A7s = new RooRealVar("A7s","A_{7}^{S}",0,-1,1); + RooRealVar* A8s = new RooRealVar("A8s","A_{8}^{S}",0,-1,1); + + // only used if Z is considered + RooRealVar* AZ = new RooRealVar("AZ","A_{Z}",0.05,0.,1.); + + RooCategory sample ("sample", "sample"); + for (unsigned int iy = 0; iy < years.size(); iy++) { + year.clear(); year.assign(Form("%i",years[iy])); + all_years += year; + + isample.clear(); isample.assign( Form("%i",nSample) ); + sample.defineType(("data"+year+"_subs"+isample).c_str()); + + } + + // Construct a simultaneous pdf using category sample as index + RooSimultaneous* simPdf = new RooSimultaneous("simPdf", "simultaneous pdf", sample); + RooSimultaneous* simPdf_penalty = new RooSimultaneous("simPdf_penalty", "simultaneous pdf with penalty term", sample); + + // Define boundary check (returning 0 in physical region and 1 outside) + BoundCheck* boundary = new BoundCheck("bound","Physical region",*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p); + + // Define boundary distance calculator + BoundDist* bound_dist = new BoundDist("bound","Physical region",*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p,true,0,false); + + // Define penalty term (parameters set to zero and will be set sample-by-sample) + Penalty* penTerm = new Penalty("penTerm","Penalty term",*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p,0,0,0,0); + + // Random generators + RooRandom::randomGenerator()->SetSeed(1); + + // + // retrieve Z4430 model workspace from file + // only a single year sample of Z is available, therefore the same parametrisation is used for the 3 years of data taking + // + RooAbsPdf* Z4430_ang_pdf = 0; + RooArgSet* Z4430_ang_params = 0; + RooAbsPdf* Z4430_mass_pdf = 0; + RooArgSet* Z4430_mass_params = 0; + if (q2Bin==6 && fitOption>0){ + string filename_Z4430 = "HistZ4430.root"; + if(!localFiles) filename_Z4430 = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/Zmodel/"+filename_Z4430; + if (!(retrieveWorkspace(filename_Z4430, wsp_Z4430, "wZ4430"))) return; + + // retrieve angular component of the Z pdf + Z4430_ang_pdf = (RooAbsPdf*) wsp_Z4430[0]->pdf("Z4430_ang_pdf"); + if(!Z4430_ang_pdf){ + std::cout << "Z4430_ang_pdf not found!!!\n" << std::endl; + wsp_Z4430[0]->allPdfs().Print(); + exit(1); + } + Z4430_ang_params = (RooArgSet*) Z4430_ang_pdf->getParameters(RooArgSet(*ctK,*ctL,*phi)); + auto iter_z = Z4430_ang_params->createIterator(); + RooRealVar* ivar_z = (RooRealVar*)iter_z->Next(); + int iii=1; + while (ivar_z) { + ivar_z->setConstant(true); + cout<GetName(),ivar_z->getVal())<Next(); + iii++; + } + // retrieve mass component of the Z pdf + Z4430_mass_pdf = (RooAbsPdf*) wsp_Z4430[0]->pdf("Z4430_mass_pdf"); + if (!Z4430_mass_pdf){ + std::cout << "Z4430_mass_pdf not found!!!\n" << std::endl; + wsp_Z4430[0]->allPdfs().Print(); + exit(1); + } + Z4430_mass_params = (RooArgSet*) Z4430_mass_pdf->getParameters(RooArgSet(*mass)); + auto iter_z_mass = Z4430_mass_params->createIterator(); + RooRealVar* ivar_z_mass = (RooRealVar*)iter_z_mass->Next(); + iii=1; + while (ivar_z_mass) { + ivar_z_mass->setConstant(true); + if (iii==3 && fitOption==2) ivar_z_mass->setConstant(false); + cout << Form("mass par %d) %s = %f\n", iii, ivar_z_mass->GetName(), ivar_z_mass->getVal()) << endl; + ivar_z_mass = (RooRealVar*) iter_z_mass->Next(); + iii++; + } + } + + // loop on the various datasets + for (unsigned int iy = 0; iy < years.size(); iy++) { + year.clear(); year.assign(Form("%i",years[iy])); + string filename_data = Form("recoTOYDataset_b%i_%i.root", q2Bin, years[iy]); + if (!localFiles) filename_data = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/toy-datasets/" + filename_data; + + // import data (or MC as data proxy) + if (!(retrieveWorkspace(filename_data, wsp, Form("ws_b%ip0", q2Bin )))) return; + + // import KDE efficiency histograms and partial integral histograms + string filename = Form((parity==0 ? "KDEeff_b%i_ev_%i.root" : "KDEeff_b%i_od_%i.root"),q2Bin,years[iy]); + if (!localFiles) { + if (XGBv<1) filename = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/eff/" + filename; + else filename = Form("/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/eff-XGBv%i/",XGBv) + filename; + } + fin_eff.push_back( new TFile( filename.c_str(), "READ" )); + if ( !fin_eff[iy] || !fin_eff[iy]->IsOpen() ) { + cout<<"File not found: "<Get(effCString.c_str())); + effWHist.push_back( (TH3D*)fin_eff[iy]->Get(effWString.c_str())); + if ( !effCHist[iy] || effCHist[iy]->IsZombie() || !effWHist[iy] || effWHist[iy]->IsZombie() ) { + cout<<"Efficiency histogram "<< effCString <<" or " << effWString << " not found in file: "<< filename <Get(intCHistString.c_str())); + intWHist.push_back( (TH1D*)fin_eff[iy]->Get(intWHistString.c_str())); + intCVec.push_back (vector (0)); + intWVec.push_back (vector (0)); + if ( !intCHist[iy] || intCHist[iy]->IsZombie() || !intWHist[iy] || intWHist[iy]->IsZombie() ) { + cout << "Integral histogram " << intCHistString <<" or " << intWHistString << " not found in file: "<< filename << endl << "Abort" << endl; + return; + } else if ( strcmp( intCHist[iy]->GetTitle(), effCHist[iy]->GetTitle() ) || strcmp( intWHist[iy]->GetTitle(), effWHist[iy]->GetTitle() )) { + // if the eff_config tag is different between efficiency and precomputed-integral means that they are inconsistent + cout << "Integral histograms are incoherent with efficiency in file: " << filename << endl; + cout << "Efficiency (CT) conf: " << effCHist[iy]->GetTitle() <GetTitle() <GetTitle() <GetTitle() <GetNbinsX(); ++i) { + intCVec[iy].push_back(intCHist[iy]->GetBinContent(i)); + } + for (int i=1; i<=intWHist[iy]->GetNbinsX(); ++i) { + intWVec[iy].push_back(intWHist[iy]->GetBinContent(i)); + } + } + + + // create roodataset (in case data-like option is selected, only import the correct % of data) + if (nSample==0) + data.push_back( createDatasetInData( wsp[iy], q2Bin, observables, shortString )); + else + data.push_back( createDatasetInToy( nSample-1, 500, wsp[iy], + q2Bin, observables, shortString )); + + if (data.back().size()==0) { + cout<<"Error producing the dataset"<0 ? Form("_MCw_xgbv%i",XGBv):""); + if (!retrieveWorkspace( filename_mc_mass, wsp_mcmass, "w")) { + cout<<"Workspace not found in mass-fit file: "<loadSnapshot(Form("reference_fit_RT_%i",q2Bin)); + RooRealVar* mean_rt = new RooRealVar (Form("mean_{RT}^{%i}",years[iy]) , "massrt" , wsp_mcmass[iy]->var(Form("mean_{RT}^{%i}",q2Bin))->getVal() , 5, 6, "GeV"); + RooRealVar* sigma_rt = new RooRealVar (Form("#sigma_{RT1}^{%i}",years[iy] ), "sigmart1" , wsp_mcmass[iy]->var(Form("#sigma_{RT1}^{%i}",q2Bin))->getVal() , 0, 1, "GeV"); + RooRealVar* alpha_rt1 = new RooRealVar (Form("#alpha_{RT1}^{%i}",years[iy] ), "alphart1" , wsp_mcmass[iy]->var(Form("#alpha_{RT1}^{%i}", q2Bin))->getVal() , 0, 10 ); + RooRealVar* n_rt1 = new RooRealVar (Form("n_{RT1}^{%i}",years[iy]) , "nrt1" , wsp_mcmass[iy]->var(Form("n_{RT1}^{%i}", q2Bin))->getVal() , 0.01, 100.); + + RooAbsPdf* dcb_rt; + RooRealVar* alpha_rt2 = new RooRealVar (Form("#alpha_{RT2}^{%i}",years[iy] ), "alphart2" , 0, -10, 10 ); + RooRealVar* n_rt2 = new RooRealVar (Form("n_{RT2}^{%i}",years[iy]) , "nrt2" , 0.01, 0.01, 100.); + RooRealVar* sigma_rt2 = new RooRealVar (Form("#sigma_{RT2}^{%i}",years[iy] ), "sigmaRT2" , 0 , 0, 0.12, "GeV"); + RooRealVar* f1rt = new RooRealVar (Form("f^{RT%i}",years[iy]) , "f1rt" , 0 , 0., 1.); + if (q2Bin != 7){ + alpha_rt2 -> setVal(wsp_mcmass[iy]->var(Form("#alpha_{RT2}^{%i}", q2Bin))->getVal() ); + n_rt2 -> setVal(wsp_mcmass[iy]->var(Form("n_{RT2}^{%i}", q2Bin) )->getVal() ); + } + + if (q2Bin >= 4){ + sigma_rt2-> setVal(wsp_mcmass[iy]->var(Form("#sigma_{RT2}^{%i}",q2Bin))->getVal() ); + f1rt -> setVal(wsp_mcmass[iy]->var(Form("f^{RT%i}", q2Bin))->getVal() ); + if (q2Bin < 7) { + if (nSample==0) dcb_rt = createRTMassShape(q2Bin, mass, mean_rt, sigma_rt, sigma_rt2, alpha_rt1, alpha_rt2, n_rt1, n_rt2 ,f1rt, wsp_mcmass[iy], years[iy], true, c_vars_rt, c_pdfs_rt ); + else dcb_rt = createRTMassShape(q2Bin, mass, mean_rt, sigma_rt, sigma_rt2, alpha_rt1, alpha_rt2, n_rt1, n_rt2 ,f1rt, wsp_mcmass[iy], years[iy], true, c_vars_rt, c_pdfs_rt, q2stat ); + } else { + if (nSample==0) dcb_rt = createRTMassShape(q2Bin, mass, mean_rt, sigma_rt, sigma_rt2, alpha_rt1, n_rt1 ,f1rt, q2Bin, wsp_mcmass[iy], years[iy], true, c_vars_rt, c_pdfs_rt ); + else dcb_rt = createRTMassShape(q2Bin, mass, mean_rt, sigma_rt, sigma_rt2, alpha_rt1, n_rt1 ,f1rt, q2Bin, wsp_mcmass[iy], years[iy], true, c_vars_rt, c_pdfs_rt, q2stat ); + } + } + else{ + alpha_rt2->setRange(0,10); + if (nSample==0) dcb_rt = createRTMassShape(q2Bin, mass, mean_rt, sigma_rt, alpha_rt1, alpha_rt2, n_rt1, n_rt2 , wsp_mcmass[iy], years[iy], true, c_vars_rt, c_pdfs_rt ); + else dcb_rt = createRTMassShape(q2Bin, mass, mean_rt, sigma_rt, alpha_rt1, alpha_rt2, n_rt1, n_rt2 , wsp_mcmass[iy], years[iy], true, c_vars_rt, c_pdfs_rt, q2stat ); + } + + // create constrained PDF for RT mass + RooArgList constr_rt_list = RooArgList(c_pdfs_rt); + constr_rt_list.add(*dcb_rt); + RooProdPdf * c_dcb_rt = new RooProdPdf(("c_dcb_rt_"+year).c_str(), ("c_dcb_rt_"+year).c_str(), constr_rt_list ); + c_vars.add(c_vars_rt); + + /// create WT component + wsp_mcmass[iy]->loadSnapshot(Form("reference_fit_WT_%i",q2Bin)); + + RooRealVar* sigma_wt = new RooRealVar (Form("#sigma_{WT1}^{%i}",years[iy]) , "sigmawt" , wsp_mcmass[iy]->var(Form("#sigma_{WT1}^{%i}", q2Bin))->getVal() , 0, 1, "GeV"); + RooRealVar* alpha_wt1 = new RooRealVar (Form("#alpha_{WT1}^{%i}",years[iy] ) , "alphawt1" , wsp_mcmass[iy]->var(Form("#alpha_{WT1}^{%i}", q2Bin))->getVal() , 0, 10 ); + RooRealVar* alpha_wt2 = new RooRealVar (Form("#alpha_{WT2}^{%i}",years[iy] ) , "alphawt2" , wsp_mcmass[iy]->var(Form("#alpha_{WT2}^{%i}", q2Bin))->getVal() , 0, 10 ); + RooRealVar* n_wt1 = new RooRealVar (Form("n_{WT1}^{%i}",years[iy]) , "nwt1" , wsp_mcmass[iy]->var(Form("n_{WT1}^{%i}", q2Bin))->getVal() , 0.01, 100.); + RooRealVar* n_wt2 = new RooRealVar (Form("n_{WT2}^{%i}",years[iy]) , "nwt2" , wsp_mcmass[iy]->var(Form("n_{WT2}^{%i}", q2Bin))->getVal() , 0.01, 100.); + + double mean_wt_val = wsp_mcmass[iy]->var(Form("mean_{WT}^{%i}", q2Bin))->getVal(); + double mean_wt_err = wsp_mcmass[iy]->var(Form("mean_{WT}^{%i}", q2Bin))->getError(); + + RooAbsPdf* dcb_wt; + double deltaPeakValue = mean_rt->getVal()-mean_wt_val; + double mean_rt_err = wsp_mcmass[iy]->var(Form("mean_{RT}^{%i}", q2Bin))->getError(); + double deltaPeakError = sqrt(mean_rt_err*mean_rt_err+mean_wt_err*mean_wt_err); + + RooRealVar* deltaPeakVar = new RooRealVar ( Form("deltaPeakVar^{%i}", years[iy]),Form("deltaPeakVar^{%i}", years[iy]), deltaPeakValue, 0., 0.2) ; + RooGaussian* c_deltaPeaks = new RooGaussian(Form("deltaPeaks^{%i}" , years[iy]) , "c_deltaPeaks", *deltaPeakVar, RooConst( deltaPeakValue ), RooConst(deltaPeakError )); // value to be checked + RooFormulaVar*mWT_data = new RooFormulaVar(Form("mWT_data^{%i}",years[iy]), "@0 + @1", RooArgList(*mean_rt, *deltaPeakVar)); + RooRealVar* mean_wt = (RooRealVar*) mWT_data; + cout << "deltaPeak constraint built, of value " << deltaPeakValue << " +/- " << deltaPeakError << endl; + + c_pdfs_wt.add(*c_deltaPeaks); + c_vars_wt.add(*deltaPeakVar); + + if (nSample==0) dcb_wt = createWTMassShape(q2Bin, mass, mean_wt, sigma_wt, alpha_wt1, alpha_wt2, n_wt1, n_wt2 , wsp_mcmass[iy], years[iy], true, c_vars_wt, c_pdfs_wt ); + else dcb_wt = createWTMassShape(q2Bin, mass, mean_wt, sigma_wt, alpha_wt1, alpha_wt2, n_wt1, n_wt2 , wsp_mcmass[iy], years[iy], true, c_vars_wt, c_pdfs_wt, q2stat ); + + // create constrained PDF for WT mass + RooArgList constr_wt_list = RooArgList(c_pdfs_wt); + constr_wt_list.add(*dcb_wt); + RooProdPdf * c_dcb_wt = new RooProdPdf(("c_dcb_wt_"+year).c_str(), ("c_dcb_wt_"+year).c_str(), constr_wt_list ); + c_vars.add(c_vars_wt); + + // As the signal PDF is written as [ CT + mFrac * MT ] (see the PdfSigAngMass class), + // the mFrac parameter represents the fraction between mistagged and correctly-tagged events (n_MT/n_CT) + // Also, since the integral of the efficiencies contains the information on the mistag fraction in MC + // this parameter represents the ratio between the fitted mFrac and the MC-based one ( n_MT_data / n_CT_data * n_CT_MC / n_MT_MC ) + RooRealVar* mFrac = new RooRealVar(Form("f_{M}^{%i}",years[iy]),"mistag fraction",1, 0, 15); + + // The values of fM_sigmas are computed on data-like MC subsamples as the fluctuation of the fraction of mistagged events ( n_MT/(n_MT+n_CT) ) + // this fluctuation needs to be propagated on the quantity of mFrac, to define a Gaussian contraint + double nrt_mc = wsp_mcmass[iy]->var(Form("nRT_%i",q2Bin))->getVal(); + double nwt_mc = wsp_mcmass[iy]->var(Form("nWT_%i",q2Bin))->getVal(); + double fraction = nwt_mc / (nrt_mc + nwt_mc); + // Propagation: sigma(mFrac) = sigma(n_MT/n_CT) * n_CT/n_MT + // = sigma(fraction)/(1-fraction)^2 * (1-fraction)/fraction + // = sigma(fraction) / fraction / (1-fraction) + double frac_sigma = fM_sigmas[years[iy]][q2Bin]/fraction/(1-fraction); + if (nSample!=0) frac_sigma = fM_sigmas[years[iy]][q2stat]/fraction/(1-fraction); + + double c_fm_val = 1; + if (evalMistagSys){ + if (iy==0) c_fm_val = 1.07; + else if (iy==1) c_fm_val = 1.02; + else if (iy==2) c_fm_val = 1.03; + } + RooGaussian* c_fm = new RooGaussian(Form("c_fm^{%i}",years[iy]) , "c_fm" , *mFrac, + RooConst(c_fm_val) , + RooConst(frac_sigma) + ); + + cout << "mFrac = " << fraction << " +/- " << fM_sigmas[years[iy]][q2Bin] << " ---> R = 1 +/- " << frac_sigma << endl; + c_vars.add(*mFrac); + + // Angular Component + RooAbsReal* ang_rt = new ShapeSigAng(("PDF_sig_ang_rt_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_rt_"+year).c_str(), + *ctK,*ctL,*phi, + *Fl,*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p, + *Fs,*As,*A4s,*A5s,*A7s,*A8s, + *effC[iy], intCVec[iy], + true + ); + + RooAbsReal* ang_wt = new ShapeSigAng(("PDF_sig_ang_wt_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_wt_"+year).c_str(), + *ctK,*ctL,*phi, + *Fl,*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p, + *Fs,*As,*A4s,*A5s,*A7s,*A8s, + *effW[iy], intWVec[iy], + false + ); + + // Angular * mass component + PdfSigAngMass* pdf_sig_ang_mass = nullptr; + PdfSigAngMass* pdf_sig_ang_mass_penalty = nullptr; + if (q2Bin < 5) { + pdf_sig_ang_mass = new PdfSigAngMass( Form(sigpdfname.c_str(),years[iy]), + Form(sigpdfname.c_str(),years[iy]), + *ctK,*ctL,*phi,*mass, + *mean_rt, *sigma_rt, *alpha_rt1, *alpha_rt2, *n_rt1, *n_rt2, + *mean_wt, *sigma_wt, *alpha_wt1, *alpha_wt2, *n_wt1, *n_wt2, + *mFrac, *c_fm, + *ang_rt, *ang_wt, + *c_dcb_rt, *c_dcb_wt + ); + + pdf_sig_ang_mass_penalty = new PdfSigAngMass( Form((sigpdfname+"_penalty").c_str(),years[iy]), + Form((sigpdfname+"_penalty").c_str(),years[iy]), + *ctK,*ctL,*phi,*mass, + *mean_rt, *sigma_rt, *alpha_rt1, *alpha_rt2, *n_rt1, *n_rt2, + *mean_wt, *sigma_wt, *alpha_wt1, *alpha_wt2, *n_wt1, *n_wt2, + *mFrac, *c_fm, + *penTerm, + *ang_rt, *ang_wt, + *c_dcb_rt, *c_dcb_wt + ); + } + else if (q2Bin==7){ + pdf_sig_ang_mass = new PdfSigAngMass( Form(sigpdfname.c_str(),years[iy]), + Form(sigpdfname.c_str(),years[iy]), + *ctK,*ctL,*phi,*mass, + *mean_rt, *sigma_rt, *sigma_rt2, *alpha_rt1, *n_rt1, *f1rt, + *mean_wt, *sigma_wt, *alpha_wt1, *alpha_wt2, *n_wt1, *n_wt2, + *mFrac, *c_fm, + *ang_rt, *ang_wt, + *c_dcb_rt, *c_dcb_wt + ); + + pdf_sig_ang_mass_penalty = new PdfSigAngMass( Form((sigpdfname+"_penalty").c_str(),years[iy]), + Form((sigpdfname+"_penalty").c_str(),years[iy]), + *ctK,*ctL,*phi,*mass, + *mean_rt, *sigma_rt, *sigma_rt2, *alpha_rt1, *n_rt1, *f1rt, + *mean_wt, *sigma_wt, *alpha_wt1, *alpha_wt2, *n_wt1, *n_wt2, + *mFrac, *c_fm, + *penTerm, + *ang_rt, *ang_wt, + *c_dcb_rt, *c_dcb_wt + ); + } + else { + pdf_sig_ang_mass = new PdfSigAngMass( Form(sigpdfname.c_str(),years[iy]), + Form(sigpdfname.c_str(),years[iy]), + *ctK,*ctL,*phi,*mass, + *mean_rt, *sigma_rt, *sigma_rt2, *alpha_rt1, *alpha_rt2, *n_rt1, *n_rt2, *f1rt, + *mean_wt, *sigma_wt, *alpha_wt1, *alpha_wt2, *n_wt1, *n_wt2, + *mFrac, *c_fm, + *ang_rt, *ang_wt, + *c_dcb_rt, *c_dcb_wt + ); + + pdf_sig_ang_mass_penalty = new PdfSigAngMass( Form((sigpdfname+"_penalty").c_str(),years[iy]), + Form((sigpdfname+"_penalty").c_str(),years[iy]), + *ctK,*ctL,*phi,*mass, + *mean_rt, *sigma_rt, *sigma_rt2, *alpha_rt1, *alpha_rt2, *n_rt1, *n_rt2, *f1rt, + *mean_wt, *sigma_wt, *alpha_wt1, *alpha_wt2, *n_wt1, *n_wt2, + *mFrac, *c_fm, + *penTerm, + *ang_rt, *ang_wt, + *c_dcb_rt, *c_dcb_wt + ); + } + + // Add constraint on mFrac to the pdf + auto pdf_sig_ang_mass_mfc = new RooProdPdf(("PDF_sig_ang_mass_mfc_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_mass_mfc_"+year).c_str(), + *pdf_sig_ang_mass, + *c_fm); + auto pdf_sig_ang_mass_penalty_mfc = new RooProdPdf(("PDF_sig_ang_mass_penalty_mfc_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_mass_penalty_mfc_"+year).c_str(), + *pdf_sig_ang_mass_penalty, + *c_fm); + + + //// Background components //// + // Read angular pdf for sidebands from external file + string filename_sb = Form("savesb_%i_b%i.root", years[iy], q2Bin); + if (!localFiles) filename_sb = "/eos/user/a/aboletti/BdToKstarMuMu/fileIndex/sidebands/" + filename_sb; + if (!(retrieveWorkspace(filename_sb, wsp_sb, "wsb"))) return; + + RooBernsteinSideband* bkg_ang_pdf = (RooBernsteinSideband*) wsp_sb[iy]->pdf(Form("BernSideBand_bin%i_%i", q2Bin, years[iy])); + RooArgSet* bkg_ang_params = (RooArgSet*)bkg_ang_pdf->getParameters(observables); + auto iter = bkg_ang_params->createIterator(); + RooRealVar* ivar = (RooRealVar*)iter->Next(); + while (ivar) { + ivar->setConstant(true); + ivar = (RooRealVar*) iter->Next(); + } + + // create exponential mass pdf for background + // in Jpsi bin, if fitOption > 0 -> bernstein polynbomial + RooRealVar* slope = new RooRealVar(Form("slope^{%i}",years[iy]), Form("slope^{%i}",years[iy]) , -5., -10., 0.); + RooAbsPdf* bkg_mass_pdf = 0; + RooExponential* bkg_exp_pdf = 0; + double pol_bmax =1.; + RooRealVar* b0_bkg_mass = new RooRealVar(Form("b0_bkg_mass-%i",years[iy]) , Form("b0_bkg_mass-%i",years[iy]) , pol_bmax ); + RooRealVar* b1_bkg_mass = new RooRealVar(Form("b1_bkg_mass-%i",years[iy]) , Form("b1_bkg_mass-%i",years[iy]) , 0.1, 0., pol_bmax); + RooRealVar* b2_bkg_mass = new RooRealVar(Form("b2_bkg_mass-%i",years[iy]) , Form("b2_bkg_mass-%i",years[iy]) , 0.1, 0., pol_bmax); + RooRealVar* b3_bkg_mass = new RooRealVar(Form("b3_bkg_mass-%i",years[iy]) , Form("b3_bkg_mass-%i",years[iy]) , 0.0 ); + RooRealVar* b4_bkg_mass = new RooRealVar(Form("b4_bkg_mass-%i",years[iy]) , Form("b4_bkg_mass-%i",years[iy]) , 0.1 , 0., pol_bmax); + if (q2Bin==4 && fitOption>0){ + b0_bkg_mass->setConstant(true); + b3_bkg_mass->setConstant(true); + bkg_mass_pdf = new RooBernstein(Form("bkg_mass1_%i",years[iy]),Form("bkg_mass1_%i",years[iy]), *mass, RooArgList(*b0_bkg_mass, *b1_bkg_mass, *b2_bkg_mass, *b3_bkg_mass,* b4_bkg_mass)); + } else { + bkg_exp_pdf = new RooExponential(Form("bkg_mass1_%i",years[iy]), Form("bkg_mass1_%i",years[iy]) , *mass, *slope ); + if (q2Bin==3||q2Bin==5) { + bkg_mass_pdf = createPdfCuts(q2Bin,years[iy], mass, slope); + }else{ + bkg_mass_pdf = bkg_exp_pdf; + } + + } + + RooProdPdf* bkg_pdf = new RooProdPdf(Form(bkgpdfname.c_str(),years[iy]), + Form(bkgpdfname.c_str(),years[iy]), + RooArgList(*bkg_ang_pdf,*bkg_mass_pdf)); + + + //// Finally build full pdf, including or not Z component //// + RooAddPdf* full_pdf = 0; + RooAddPdf* full_pdf_penalty = 0; + + // fraction of signal and bkg pdfs + RooRealVar *fsig = new RooRealVar( ("fsig_"+shortString+"_"+year).c_str(), ("fsig_"+shortString+"_"+year).c_str(),0,1 ); + + + // special case of bin 6, include Z component in the signal model + if(q2Bin==6 && fitOption>0){ + RooProdPdf* Z4430_pdf = 0; + if(fitOption==3){ + std::cout << Form("Warning! Z(4430) Mass model from Signal MC Fits for Year=%i", years[iy]) << std::endl; + RooConstVar* FracZ4430WT = new RooConstVar(Form("FracZ4430WT^{%i}",years[iy]),"FracZ4430WT", fraction); + RooProduct* mass_Frac = new RooProduct( Form("mass_Frac_%i", years[iy]), Form("mass_Frac_%i", years[iy]), RooArgList(*FracZ4430WT, *mFrac)); + RooAddPdf* Mass_All = new RooAddPdf( Form("Mass_All_%i", years[iy]), Form("Mass_All_%i", years[iy]), RooArgList(*c_dcb_wt,*c_dcb_rt), *mass_Frac); + Z4430_pdf = new RooProdPdf( Form("Z4430_pdf_%i", years[iy]), Form("Z4430_pdf_%i", years[iy]), RooArgList(*Z4430_ang_pdf, *Mass_All, *c_fm) ); + }else{ + std::cout << Form("Warning! Z(4430) Mass model from Z(4430) MC Fit for Year=%i", years[iy]) << std::endl; + Z4430_pdf = new RooProdPdf( Form("Z4430_pdf_%i",years[iy]), Form("Z4430_pdf_%i", years[iy]), RooArgList(*Z4430_ang_pdf, *Z4430_mass_pdf) ); + } + + RooAddPdf* pdf_z_sig_ang_mass_mfc = new RooAddPdf(("PDF1_sig_ang_mass_"+shortString+"_"+year).c_str(), + ("PDF1_sig_ang_mass_"+year).c_str(), + RooArgList(*Z4430_pdf,*pdf_sig_ang_mass_mfc), + RooArgList( *AZ) + ); + + RooAddPdf* pdf_z_sig_ang_mass_mfc_penalty = new RooAddPdf(("PDF1_sig_ang_mass_penalty _"+shortString+"_"+year).c_str(), + ("PDF1_sig_ang_mass_penalty _"+year).c_str(), + RooArgList(*Z4430_pdf,*pdf_sig_ang_mass_penalty_mfc), + RooArgList( *AZ) + ); + + full_pdf = new RooAddPdf( ("PDF_sig_ang_fullAngularMass_bkg_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_fullAngularMass_bkg_"+shortString+"_"+year).c_str(), + RooArgList(*pdf_z_sig_ang_mass_mfc, *bkg_pdf), + RooArgList(*fsig) + ); + full_pdf_penalty = new RooAddPdf(("PDF_sig_ang_fullAngularMass_bkg_penalty_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_fullAngularMass_bkg_penalty_"+shortString+"_"+year).c_str(), + RooArgList(*pdf_z_sig_ang_mass_mfc_penalty, *bkg_pdf), + RooArgList(*fsig) + ); + + } + else { + + full_pdf = new RooAddPdf(("PDF_sig_ang_fullAngularMass_bkg_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_fullAngularMass_bkg_"+shortString+"_"+year).c_str(), + RooArgList(*pdf_sig_ang_mass_mfc, *bkg_pdf), + RooArgList(*fsig) + ); + + full_pdf_penalty = new RooAddPdf(("PDF_sig_ang_fullAngularMass_bkg_penalty_"+shortString+"_"+year).c_str(), + ("PDF_sig_ang_fullAngularMass_bkg_penalty_"+shortString+"_"+year).c_str(), + RooArgList(*pdf_sig_ang_mass_penalty_mfc, *bkg_pdf), + RooArgList(*fsig) + ); + } + PDF_sig_ang_mass_bkg.push_back(full_pdf); + PDF_sig_ang_mass_bkg_penalty.push_back(full_pdf_penalty); + + // insert sample in the category map, to be imported in the combined dataset + // and associate model with the data + if ( !data[iy][0] || data[iy][0]->IsZombie() ) { + cout<<"Dataset " << nSample << " not found in file: "<(("data"+year+Form("_subs%d",nSample)).c_str(), data[iy][0]) ); + simPdf ->addPdf(*PDF_sig_ang_mass_bkg[iy], ("data"+year+Form("_subs%d",nSample)).c_str()); + simPdf_penalty->addPdf(*PDF_sig_ang_mass_bkg_penalty[iy], ("data"+year+Form("_subs%d",nSample)).c_str()); + + } + + TFile* fout = 0; + std::string addInfo = evalMistagSys ? "_mistagSyst_exp" : ""; + + if (save>0) fout = new TFile(("simFitResults4d/simFitResult_toy_fullAngularMass_Swave_" + all_years + stat + Form("_b%i%s_unbl%i%s.root", q2Bin, XGBstr.c_str(), unblind, addInfo.c_str())).c_str(),"RECREATE"); + RooWorkspace* wsp_out = 0; + + // save initial par values + RooArgSet *params = (RooArgSet *)simPdf->getParameters(observables); + RooArgSet* savedParams = (RooArgSet *)params->snapshot() ; + // Construct combined dataset in (x,sample) + RooDataSet allcombData ("allcombData", "combined data", + reco_vars, + Index(sample), + Import(map)); + RooDataSet* combData = 0; + + // Results' containers + double fitTime, imprTime, minTime; + double co1, co4, co5; + double boundDistFit, boundDist; + bool boundCheck, convCheck, usedPenalty; + + RooArgList pars (*Fl,*P1,*P2,*P3,*P4p,*P5p,*P6p,*P8p); + RooArgList sPars (*Fs,*As,*A4s,*A5s,*A7s,*A8s); + + // TTree with the MINOS output + vector vResult (pars.getSize()); + vector vConfInterLow (pars.getSize()); + vector vConfInterHigh (pars.getSize()); + vector vFreeResult (pars.getSize()); + vector vFreeConfInterLow (pars.getSize()); + vector vFreeConfInterHigh (pars.getSize()); + if (save>0) fout->cd(); + TTree* fitResultsTree = new TTree("fitResultsTree","fitResultsTree"); + for (int iPar = 0; iPar < pars.getSize(); ++iPar) { + RooRealVar* par = (RooRealVar*)pars.at(iPar); + fitResultsTree->Branch(Form("%s_low",par->GetName()),&vConfInterLow[iPar]); + fitResultsTree->Branch(Form("%s_high",par->GetName()),&vConfInterHigh[iPar]); + fitResultsTree->Branch(Form("%s_best",par->GetName()),&vResult[iPar]); + fitResultsTree->Branch(Form("%s_free_low",par->GetName()),&vFreeConfInterLow[iPar]); + fitResultsTree->Branch(Form("%s_free_high",par->GetName()),&vFreeConfInterHigh[iPar]); + fitResultsTree->Branch(Form("%s_free_best",par->GetName()),&vFreeResult[iPar]); + } + fitResultsTree->Branch("fitTime",&fitTime); + fitResultsTree->Branch("imprTime",&imprTime); + fitResultsTree->Branch("minTime",&minTime); + fitResultsTree->Branch("co1",&co1); + fitResultsTree->Branch("co4",&co4); + fitResultsTree->Branch("co5",&co5); + fitResultsTree->Branch("boundDist",&boundDist); + fitResultsTree->Branch("boundDistFit",&boundDistFit); + fitResultsTree->Branch("boundCheck",&boundCheck); + fitResultsTree->Branch("convCheck",&convCheck); + fitResultsTree->Branch("usedPenalty",&usedPenalty); + + // Timer for fitting time + TStopwatch subTime; + + // counters to monitor results' status + int cnt[9]; + for (int iCnt=0; iCnt<9; ++iCnt) cnt[iCnt] = 0; + + Fitter* fitter = 0; + vector vFitter (0); + + { + uint is = nSample; + + string samplename = Form("_subs%d", is); + samplename = "data%i" + samplename; + string the_cut = Form(("sample==sample::"+samplename).c_str(),years[0]); + if (years.size() > 1){ + for (unsigned int iy=1; iy < years.size(); iy++){ + the_cut = the_cut + Form("|| sample==sample::data%d_subs%d", years[iy], is); + } + } + + combData = (RooDataSet*)allcombData.reduce(Cut(the_cut.c_str())); + if (nSample>0) cout<<"Fitting data subsample "<numEntries()<<" entries"<numEntries()<<" entries"<numEntries(); + penTerm->setPower(power/combEntries); + + // to start the fit, parameters are restored to the center of the parameter space + *params = *savedParams ; + + // run the fit + fitter = new Fitter (Form("fitter%i",is),Form("fitter%i",is),pars,sPars,combData,simPdf,simPdf_penalty,boundary,bound_dist,penTerm,&c_vars); + vFitter.push_back(fitter); + + fitter->setUnbl(unblind); + + bool runPostFitSteps = true; + if (nSample==0 && (q2Bin==4 || q2Bin==6)) { + // define if run improvAng and minosAng + runPostFitSteps = false; + fitter->setNCPU(8); + // for bin 4, do not run the penalty even if needed + if (q2Bin==4) fitter->runSimpleFit = true; + } + +// if (q2Bin==3 || q2Bin==5){ +// fitter->DisableConstOpt(); +// std::cout<fit(); + subTime.Stop(); + + fitTime=subTime.CpuTime(); + cout<<(fitter->runSimpleFit?"Fit time: ":"Fit+boundDist time: ")<getValV() == 0; + + if (unblind>3) fitter->result()->Print("v"); + + if (fitter->runSimpleFit) boundDistFit = boundDist = -1; + else boundDistFit = boundDist = fitter->boundDist; + + usedPenalty = fitter->usedPenalty; + + if (usedPenalty) { + // save coefficient values + co1 = fitter->coeff1; + co4 = fitter->coeff4; + co5 = fitter->coeff5; + + TStopwatch improvTime; + improvTime.Start(true); + if (runPostFitSteps) fitter->improveAng(1,300000); + improvTime.Stop(); + imprTime = improvTime.CpuTime(); + cout<<"Improv time: "<boundDist; + + } + + // run MINOS error + TStopwatch minosTime; + minosTime.Start(true); + + if (runPostFitSteps) fitter->MinosAng(); + + minosTime.Stop(); + minTime = minosTime.CpuTime(); + cout<<"MINOS errors computed in "<vResult[iPar]; + vFreeResult[iPar] = fitter->vFreeFitResult[iPar]; + if (runPostFitSteps) { + vConfInterLow[iPar] = fitter->vConfInterLow[iPar] - vResult[iPar]; + vConfInterHigh[iPar] = fitter->vConfInterHigh[iPar] - vResult[iPar]; + } else { + vConfInterLow[iPar] = fitter->vFitErrLow[iPar]; + vConfInterHigh[iPar] = fitter->vFitErrHigh[iPar]; + } + vFreeConfInterLow[iPar] = fitter->vFreeFitErrLow[iPar]; + vFreeConfInterHigh[iPar] = fitter->vFreeFitErrHigh[iPar]; + } + fitResultsTree->Fill(); + + if (plot && !multiSample && unblind>2) { + + string plotString = shortString + "_" + all_years + stat + XGBstr + Form("_unbl%i",unblind); + addInfo = evalMistagSys ? "_mistagSyst_exp" : ""; + string plotname = "plotSimFit4d_d/simFitResult_toy_fullAngularMass_Swave_" + plotString + addInfo + ".pdf"; + fitter->plotSimFitProjections(plotname.c_str(),{samplename,sigpdfname,bkgpdfname},years,true); + + } + + if (save>1 && (q2Bin!=4 || years.size()<3 || nSample>0) && unblind>1) { + wsp_out = new RooWorkspace("wsp_out","wsp_out"); + if (unblind<4) + for (int iPar = 0; iPar < pars.getSize(); ++iPar) { + ((RooRealVar*)pars.at(iPar))->setVal(0); + ((RooRealVar*)pars.at(iPar))->setError(0); + } + wsp_out->import(*combData); + wsp_out->import(*simPdf); + } + + } + + // fill fit-status-dependent counters + ++cnt[8]; + int iCnt = 0; + if (!convCheck) iCnt += 4; + if (!boundCheck) iCnt += 2; + if (usedPenalty) iCnt += 1; + ++cnt[iCnt]; + + // print fit status and time + if (!boundCheck) + if (convCheck) cout<<"Converged in unphysical region"; + else cout<<"Not converged"; + else + if (convCheck) + if (usedPenalty) cout<<"Converged with penalty term with coeff: "<coeff1<<" "<coeff4<<" "<coeff5; + else cout<<"Converged without penalty"; + else cout<<"This should never be printed"; + cout<<" ("<0 && unblind>1) { + fout->cd(); + if (unblind>3) fitResultsTree->Write(); + if (wsp_out) wsp_out->Write(); + fout->Close(); + } + +} + + +void simfit_data_fullAngularMass_SwaveBin1(int q2Bin, int parity, bool multiSample, uint nSample, uint q2stat, int fitOption, int XGBv, int unblind, bool localFiles, bool plot, int save, std::vector years) +{ + if ( parity==-1 ) + for (parity=0; parity<2; ++parity) + simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, multiSample, nSample, q2stat, fitOption, XGBv, unblind, localFiles, plot, save, years); + else + simfit_data_fullAngularMass_SwaveBin(q2Bin, parity, multiSample, nSample, q2stat, fitOption, XGBv, unblind, localFiles, plot, save, years); +} + +int main(int argc, char** argv) +{ + // q2-bin format: [0-8] for one bin + // [-1] for each bin recursively + // parity format: [0] even efficiency + // [1] odd efficiency + // [-1] for each parity recursively + // unblind options: [0] no fit in signal q2 bins + // [1] fit reports only time, convergence and penalty + // [2] adds root file with RooWorkspace and blinded POI + // [3] adds fit projections plot + // [4] adds RooFitResult and MINOS errors to log, unblind POI and add TTree in root file + + int q2Bin = -1; + int parity = -1; + + if ( argc > 1 ) q2Bin = atoi(argv[1]); + if ( argc > 2 ) parity = atoi(argv[2]); + + bool multiSample = false; + uint nSample = 0; + uint q2stat = 0; + if ( argc > 3 && atoi(argv[3]) > 0 ) multiSample = true; + if ( argc > 4 ) nSample = atoi(argv[4]); + if ( argc > 5 ) q2stat = atoi(argv[5]); + + if (nSample==0) multiSample = false; + + int fitOption=0; + // possible configurations and meanings of alternative configurations: + // bin 4: fitOpt == 1 -> bkg mass pdf bernstein polynomial instead of expo + // bin 6: fitOpt == 1 -> include Z component, all pars fixed to MC + // fitOpt == 2 -> include Z component, some pars are floating + // fitOpt == 3 -> include Z component, Z mass model from the psi2S MC + if ( argc > 6 ) fitOption = atoi(argv[6]); + + int XGBv = 0; + if ( argc > 7 ) XGBv = atoi(argv[7]); + + int unblind=0; + if ( argc > 8 ) unblind = atoi(argv[8]); + + bool localFiles = false; + if ( argc > 9 && atoi(argv[9]) > 0 ) localFiles = true; + + bool plot = true; + int save = true; + + if ( argc > 10 && atoi(argv[10]) == 0 ) plot = false; + if ( argc > 11 ) save = atoi(argv[11]); + + std::vector years; + if ( argc > 12 && atoi(argv[12]) != 0 ) years.push_back(atoi(argv[12])); + else { + cout << "No specific years selected, using default: 2016" << endl; + years.push_back(2016); + } + if ( argc > 13 && atoi(argv[13]) != 0 ) years.push_back(atoi(argv[13])); + if ( argc > 14 && atoi(argv[14]) != 0 ) years.push_back(atoi(argv[14])); + + cout << "q2Bin " << q2Bin << endl; + cout << "parity " << parity << endl; + cout << "multiSample " << multiSample << endl; + cout << "nSample " << nSample << endl; + cout << "q2stat " << q2stat << endl; + cout << "fitOption " << fitOption << endl; + cout << "XGB version " << XGBv << endl; + cout << "local files " << localFiles << endl; + cout << "plot " << plot << endl; + cout << "save " << save << endl; + cout << "years[0] " << years[0] << endl; +// cout << years[1] << endl; +// cout << years[2] << endl; + + cout << "UNBLINDIG STEP " << unblind << endl; + + if ( q2Bin < -1 || q2Bin >= nBins ) return 1; + if ( parity < -1 || parity > 1 ) return 1; + + if ( q2stat < 0 || q2stat >= nBins ) return 1; + + if ( XGBv < 0 ) return 1; + + // Protectrion against accidental unblinding + if ( unblind < 1 && q2Bin != 4 && q2Bin != 6 ) { + cout<<"The analysis is blind!"< logs_simFit4d/simfit_toy_fullAngularMass_Swave_${bin}_${par}_${multi}_${nsam}_${q2stat}_${fitopt}_${XGBv}_${localFile}_${plot}_${save}_2016_2017_2018_unbl${unbl}.log & + ;; + + 1) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2016 \ + &> logs_simFit4d/simfit_toy_fullAngularMass_Swave_${bin}_${par}_${multi}_${nsam}_${q2stat}_${fitopt}_${XGBv}_${localFile}_${plot}_${save}_2016_unbl${unbl}.log & + ;; + + 2) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2017 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2017 \ + &> logs_simFit4d/simfit_toy_fullAngularMass_Swave_${bin}_${par}_${multi}_${nsam}_${q2stat}_${fitopt}_${XGBv}_${localFile}_${plot}_${save}_2017_unbl${unbl}.log & + ;; + + 3) + echo ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2018 + ./simfit_toy_fullAngularMass_Swave ${bin} ${par} ${multi} ${nsam} ${q2stat} ${fitopt} ${XGBv} ${unbl} ${localFile} ${plot} ${save} 2018 \ + &> logs_simFit4d/simfit_toy_fullAngularMass_Swave_${bin}_${par}_${multi}_${nsam}_${q2stat}_${fitopt}_${XGBv}_${localFile}_${plot}_${save}_2018_unbl${unbl}.log & + ;; + + esac + +fi diff --git a/src/BoundDist.cc b/src/BoundDist.cc index 7663d66..b6ec235 100644 --- a/src/BoundDist.cc +++ b/src/BoundDist.cc @@ -106,7 +106,7 @@ Double_t BoundDist::evalR(double Radius) const bool nextRun; - TRandom3 randGen(0); + TRandom3 randGen(1); for (int iProbe=0; iProbe(angPars.getSize(),0); vConfInterHigh = std::vector(angPars.getSize(),0); + vSResult = std::vector(sPars.getSize(),0); + vSErrLow = std::vector(sPars.getSize(),0); + vSErrHigh = std::vector(sPars.getSize(),0); + vFreeSResult = std::vector(sPars.getSize(),0); + vFreeSErrLow = std::vector(sPars.getSize(),0); + vFreeSErrHigh = std::vector(sPars.getSize(),0); + nCPU = 1; nCPU_Pen=1; unbl = 0; + printMinosScan = false; + + disableConstOpt = false; + } @@ -162,7 +235,8 @@ Int_t Fitter::fit() } RooMinimizer m(*nll) ; - m.optimizeConst (kTRUE); // do not recalculate constant terms + if (!disableConstOpt) + m.optimizeConst (kTRUE); // do not recalculate constant terms m.setOffsetting(kTRUE); // Enable internal likelihood offsetting for enhanced numeric precision. // m.setVerbose(kTRUE); m.setPrintLevel(-1); @@ -308,6 +382,7 @@ Int_t Fitter::improveAng(int seed, int nGen) double improvNLL = NLL_before; double testNLL = 0; int iImprove = 0; + std::cout << Form("NLL_before = %.5f",NLL_before) << std::endl; do { @@ -335,6 +410,7 @@ Int_t Fitter::improveAng(int seed, int nGen) ((RooRealVar*)angPars.at(iPar))->setVal(vImprovPar[iPar]); computeBoundaryDistance(); + std::cout << Form("NLL_after = %.5f",improvNLL) << std::endl; std::cout<<"Improved fit result: deltaNLL = "< "< vLastHit(0); + std::vector vLastSHit(0); for (int iPar1 = 0; iPar1 < angPars.getSize(); ++iPar1) vLastHit.push_back(vResult[iPar1]); + for (int iPar1 = 0; iPar1 < sPars.getSize(); ++iPar1) + vLastSHit.push_back(vSResult[iPar1]); TH1D* parRandomPool = 0; int nHistBins = 0; @@ -413,18 +501,35 @@ Int_t Fitter::MinosAng(int seed, int nGenMINOS) do p_test = parRandomPool->GetRandom(); while (p_test>par->getMax() || p_testgetMin()); par->setVal(p_test); + std::string nllprint = "l "; for (int iPar1 = 0; iPar1 < angPars.getSize(); ++iPar1) { - if (iPar1==iPar) continue; + if (iPar1==iPar) { + nllprint = nllprint + Form("%.4f ",p_test); + continue; + } RooRealVar* par1 = (RooRealVar*)angPars.at(iPar1); double par1val = 0; do par1val = randGenMinos.Gaus(vLastHit[iPar1],widthScale*TMath::Max(0.5*(vFitErrHigh[iPar1]-vFitErrLow[iPar1]),minParError)); while (par1val>par1->getMax() || par1valgetMin()); par1->setVal(par1val); + nllprint = nllprint + Form("%.4f ",par1val); } // check if the point is physical if (boundary->getValV()>0) continue; // get and test the local likelihood + for (int iPar1 = 0; iPar1 < sPars.getSize(); ++iPar1) { + RooRealVar* par1 = (RooRealVar*)sPars.at(iPar1); + double par1val = 0; + do par1val = randGenMinos.Gaus(vLastSHit[iPar1],widthScale*TMath::Max(0.5*(vSErrHigh[iPar1]-vSErrLow[iPar1]),minParError)); + while (par1val>par1->getMax() || par1valgetMin()); + par1->setVal(par1val); + nllprint = nllprint + Form("%.4f ",par1val); + } probedNLL = nll->getValV(); + + if (printMinosScan) + std::cout << nllprint << Form("%.3f",probedNLL) << std::endl; + if (probedNLL<=NLL_min+0.5) { p_in = p_test; if ( isErrHigh > 0 ) { if ( p_in > par->getMax()-parRandomPool->GetBinWidth(1) ) break; } @@ -433,6 +538,10 @@ Int_t Fitter::MinosAng(int seed, int nGenMINOS) RooRealVar* par1 = (RooRealVar*)angPars.at(iPar1); vLastHit[iPar1] = par1->getValV(); } + for (int iPar1 = 0; iPar1 < sPars.getSize(); ++iPar1) { + RooRealVar* par1 = (RooRealVar*)sPars.at(iPar1); + vLastSHit[iPar1] = par1->getValV(); + } if (isErrHigh>0) for (int iBin=1; iBin<=parRandomPool->FindBin(p_test); ++iBin) parRandomPool->SetBinContent(iBin,0); @@ -494,6 +603,12 @@ Int_t Fitter::MinosAng(int seed, int nGenMINOS) } + // Reset parameters to best fit result + for (int iPar = 0; iPar < angPars.getSize(); ++iPar) + ((RooRealVar*)angPars.at(iPar))->setVal(vResult[iPar]); + for (int iPar = 0; iPar < sPars.getSize(); ++iPar) + ((RooRealVar*)sPars.at(iPar))->setVal(vSResult[iPar]); + return 0; } @@ -526,6 +641,19 @@ void Fitter::fillResultContainers(bool fromImprov, bool isFree) } + if (!fromImprov) + for (int iPar = 0; iPar < sPars.getSize(); ++iPar) { + RooRealVar* par = (RooRealVar*)sPars.at(iPar); + vSResult[iPar] = par->getValV(); + vSErrLow[iPar] = par->getErrorLo(); + vSErrHigh[iPar] = par->getErrorHi(); + if (isFree) { + vFreeSResult[iPar] = vResult[iPar]; + vFreeSErrLow[iPar] = vSErrLow[iPar]; + vFreeSErrHigh[iPar] = vSErrHigh[iPar]; + } + } + } @@ -597,10 +725,9 @@ void Fitter::plotProjections(RooAbsPdf* singleYearPdf, RooAbsData* singleYearDat frames.push_back( var->frame(RooFit::Title(Form(frametitle.c_str(),var->GetTitle()))) ); singleYearData->plotOn(frames[fr], - RooFit::MarkerColor(kRed+1), - RooFit::LineColor(kRed+1), - RooFit::Binning(40), - RooFit::Name(Form("plData%i",ipad)) ); + // RooFit::MarkerColor(kRed+1), + // RooFit::LineColor(kRed+1), + RooFit::Binning(40) ); singleYearPdf->plotOn(frames[fr], RooFit::LineWidth(1), @@ -623,6 +750,10 @@ void Fitter::plotProjections(RooAbsPdf* singleYearPdf, RooAbsData* singleYearDat RooFit::Components( catnames[1].c_str() )); } + singleYearData->plotOn(frames[fr], + RooFit::Binning(40), + RooFit::Name(Form("plData%i",ipad)) ); + if (fr == 0) { leg->AddEntry(frames[fr]->findObject(Form("plData%i",ipad)), "Data", "lep"); diff --git a/src/RooBernsteinSideband.cxx b/src/RooBernsteinSideband.cxx index a3b4e02..feade9a 100644 --- a/src/RooBernsteinSideband.cxx +++ b/src/RooBernsteinSideband.cxx @@ -1,10 +1,17 @@ /***************************************************************************** - * Project: RooBernsteinSideband * + * Project: RooBernsteinSideband * * * * P.Dini fecit, Anno Domini MMXVIII * + * * * "Est modus in rebus." * * * - * Class to describe 3D angular Sidebandciency in B0->K*MuMU Analysis * + * Class to describe 3D angular Sideband in B0->K*mumu Analysis * + * * + * Ipse mutavit, Anno Domini MMXXIV * + * * + * "Et de hoc satis" * + * * + * Now with a trivial circularity condition added * * * *****************************************************************************/ @@ -52,7 +59,7 @@ ClassImp(RooBernsteinSideband); int icc=0; for(int i = 0; i <= maxDegree1 ; ++i) { for(int j = 0; j <= maxDegree2 ; ++j) { - for(int k = 0; k <= maxDegree3 ; ++k) { + for(int k = 0; k < maxDegree3 ; ++k) { printf("RooBernsteinSideband: Par(%d)=%f cosL=%d cosK=%d phi=%d\n",icc,((RooAbsReal&) _coefList[icc]).getVal(), i,j,k); icc++; @@ -141,6 +148,7 @@ Double_t RooBernsteinSideband::evaluate() const // long double sz_div = 1.0/(1.0-z); int ipar =0; + int ipa0 =0; double func =0.0; double intg_1 =0.0; @@ -164,10 +172,16 @@ Double_t RooBernsteinSideband::evaluate() const // func += ((RooAbsReal&) _coefList[ipar-1]).getVal()*sx[_maxDegree1-i]*sy[_maxDegree2-j]*sz[_maxDegree3-k]; // intg_1 +=((RooAbsReal&) _coefList[ipar-1]).getVal(); // if(fabs(((RooAbsReal&) _coefList[ipar]).getVal())>0.0){ + if (k==0) ipa0=ipar; bernknvalz = device_coeffbinomial(_maxDegree3,k)*tz*sz[_maxDegree3-k]; - func += ((RooAbsReal&) _coefList[ipar]).getVal()*bernknvalx*bernknvaly*bernknvalz; - intg_1 += ((RooAbsReal&) _coefList[ipar]).getVal(); - ipar++; + if(k==_maxDegree3){ + func += ((RooAbsReal&) _coefList[ipa0]).getVal()*bernknvalx*bernknvaly*bernknvalz; + intg_1 += ((RooAbsReal&) _coefList[ipa0]).getVal(); + }else{ + func += ((RooAbsReal&) _coefList[ipar]).getVal()*bernknvalx*bernknvaly*bernknvalz; + intg_1 += ((RooAbsReal&) _coefList[ipar]).getVal(); + ipar++; + } // } // std::cout<<"coeff p("< temp_sub_simfit_recoMC_fullAngular_oneBin.sub +nbins=$(wc -l ../confSF/KDE_SF_all.list | cut -d " " -f1) + +while read -a line; do + bin=${line[0]} + # Creation of the submit HTCondor file + cat << EOF > temp_sub_simfit_recoMC_fullAngular_oneBin${bin}.sub Executable = run_simfit_recoMC_fullAngular.sh -nsamp = ( \$(ProcId) / 6 ) + 1 -bin = \$(ProcId) % 6 -Arguments = \$INT(nsamp) \$INT(bin) +nsamp = \$(ProcId) +parity = 1 +yearConf = 0 +Arguments = \$INT(nsamp) ${bin} \$INT(parity) \$INT(yearConf) Log = logs_parSub/sub_\$(ClusterId).log -Output = logs_parSub/simfit_recoMC_fullAngular_randLikelihood_\$INT(nsamp)_\$INT(bin).out -Error = logs_parSub/simfit_recoMC_fullAngular_randLikelihood_\$INT(nsamp)_\$INT(bin).err +Output = logs_parSub/simfit_recoMC_fullAngular_\$INT(nsamp)_p\$INT(parity)_b${bin}.out +Error = logs_parSub/simfit_recoMC_fullAngular_\$INT(nsamp)_p\$INT(parity)_b${bin}.err transfer_output_files = "" +JobFlavour = "testmatch" EOF -if [ "${USER}" == "fiorendi" ]; then - echo '+AccountingGroup = "group_u_CMST3.all"'>>temp_sub_simfit_recoMC_fullAngular_oneBin.sub -fi -echo 'Queue 600'>>temp_sub_simfit_recoMC_fullAngular_oneBin.sub + if [ "${USER}" == "fiorendi" ]; then + echo '+AccountingGroup = "group_u_CMST3.all"'>>temp_sub_simfit_recoMC_fullAngular_oneBin{bin}.sub + fi + echo "Queue 100">>temp_sub_simfit_recoMC_fullAngular_oneBin${bin}.sub -# Compilation, submission and file removal -if make simfit_recoMC_fullAngular -then condor_submit temp_sub_simfit_recoMC_fullAngular_oneBin.sub -fi -rm temp_sub_simfit_recoMC_fullAngular_oneBin.sub + # Compilation, submission and file removal + if make simfit_recoMC_fullAngular + then condor_submit temp_sub_simfit_recoMC_fullAngular_oneBin${bin}.sub + fi + rm temp_sub_simfit_recoMC_fullAngular_oneBin{bin}.sub +done < ../confSF/KDE_SF_all.list diff --git a/sub_simfit_recoMC_fullAngularMass.sh b/sub_simfit_recoMC_fullAngularMass.sh index 04efaa6..3227aad 100644 --- a/sub_simfit_recoMC_fullAngularMass.sh +++ b/sub_simfit_recoMC_fullAngularMass.sh @@ -1,7 +1,7 @@ #!/bin/bash # Create directory for log files -if [ ! -d logs_parSub/xgbv8 ]; then mkdir -p logs_parSub/xgbv8; fi +if [ ! -d logs_parSub ]; then mkdir -p logs_parSub; fi nbins=$(wc -l ../confSF/KDE_SF_all.list | cut -d " " -f1) @@ -10,11 +10,13 @@ while read -a line; do # Creation of the submit HTCondor file cat << EOF > temp_sub_simfit_recoMC_fullAngularMass_oneBin${bin}.sub Executable = run_simfit_recoMC_fullAngularMass.sh +parity = 0 nsamp = \$(ProcId) -Arguments = \$INT(nsamp) ${bin} -Log = logs_parSub/xgbv8/sub_\$(ClusterId).log -Output = logs_parSub/xgbv8/simfit_recoMC_fullAngularMass_\$(ClusterId)_\$(ProcId)_\$INT(nsamp)_${bin}.out -Error = logs_parSub/xgbv8/simfit_recoMC_fullAngularMass_\$(ClusterId)_\$(ProcId)_\$INT(nsamp)_${bin}.err +yearConf = 0 +Arguments = \$INT(nsamp) ${bin} \$INT(parity) \$INT(yearConf) +Log = logs_parSub/sub_\$(ClusterId).log +Output = logs_parSub/simfit_recoMC_fullAngularMass_\$INT(nsamp)_p\$INT(parity)_b${bin}.out +Error = logs_parSub/simfit_recoMC_fullAngularMass_\$INT(nsamp)_p\$INT(parity)_b${bin}.err transfer_output_files = "" +JobFlavour = "testmatch" EOF diff --git a/sub_simfit_recoMC_fullAngularMass_toybkg.sh b/sub_simfit_recoMC_fullAngularMass_toybkg.sh index 328d6e2..902fae3 100644 --- a/sub_simfit_recoMC_fullAngularMass_toybkg.sh +++ b/sub_simfit_recoMC_fullAngularMass_toybkg.sh @@ -1,7 +1,7 @@ #!/bin/bash # Create directory for log files -if [ ! -d logs_parSub/xgbv8 ]; then mkdir -p logs_parSub/xgbv8; fi +if [ ! -d logs_parSub ]; then mkdir -p logs_parSub; fi nbins=$(wc -l ../confSF/KDE_SF_all.list | cut -d " " -f1) @@ -10,11 +10,12 @@ while read -a line; do # Creation of the submit HTCondor file cat << EOF > temp_sub_simfit_recoMC_fullAngularMass_toybkg_oneBin${bin}.sub Executable = run_simfit_recoMC_fullAngularMass_toybkg.sh +parity = 0 nsamp = \$(ProcId) -Arguments = \$INT(nsamp) ${bin} -Log = logs_parSub/xgbv8/sub_\$(ClusterId).log -Output = logs_parSub/xgbv8/simfit_recoMC_fullAngularMass_toybkg_\$INT(nsamp)_${bin}.out -Error = logs_parSub/xgbv8/simfit_recoMC_fullAngularMass_toybkg_\$INT(nsamp)_${bin}.err +Arguments = \$INT(nsamp) ${bin} \$INT(parity) +Log = logs_parSub/sub_\$(ClusterId).log +Output = logs_parSub/simfit_recoMC_fullAngularMass_toybkg_\$INT(nsamp)_p\$INT(parity)_b${bin}.out +Error = logs_parSub/simfit_recoMC_fullAngularMass_toybkg_\$INT(nsamp)_p\$INT(parity)_b${bin}.err transfer_output_files = "" +JobFlavour = "testmatch" EOF diff --git a/sub_simfit_toy_fullAngularMass_Swave.sh b/sub_simfit_toy_fullAngularMass_Swave.sh new file mode 100644 index 0000000..ea6646c --- /dev/null +++ b/sub_simfit_toy_fullAngularMass_Swave.sh @@ -0,0 +1,33 @@ +#!/bin/bash + +# Create directory for log files +if [ ! -d logs_parSub ]; then mkdir logs_parSub; fi + +# Creation of the submit HTCondor file +cat << EOF > temp_sub_simfit_toy_fullAngularMass_Swave_oneBin.sub +Executable = run_simfit_toy_fullAngularMass_Swave.sh +nsamp = \$(ProcId) + 1 +bin = 2 +Arguments = \$INT(nsamp) \$INT(bin) +Log = logs_parSub/sub_\$(ClusterId).log +Output = logs_parSub/simfit_toy_fullAngularMass_\$INT(nsamp)_\$INT(bin).out +Error = logs_parSub/simfit_toy_fullAngularMass_\$INT(nsamp)_\$INT(bin).err +transfer_output_files = "" ++JobFlavour = "testmatch" +EOF + +if [ "${USER}" == "fiorendi" ]; then + echo '+AccountingGroup = "group_u_CMST3.all"'>>temp_sub_simfit_toy_fullAngularMass_Swave_oneBin.sub +fi +echo 'Queue 500'>>temp_sub_simfit_toy_fullAngularMass_Swave_oneBin.sub + +# Compilation, submission and file removal +if make simfit_toy_fullAngularMass_Swave +then condor_submit temp_sub_simfit_toy_fullAngularMass_Swave_oneBin.sub +fi +rm temp_sub_simfit_toy_fullAngularMass_Swave_oneBin.sub + +# bin = ( \$(ProcId) / 4 ) * 2 + 4 +# yearConf = \$(ProcId) % 4 +# Output = logs_parSub/simfit_toy_fullAngularMass_Swave_\$INT(nsamp)_\$INT(bin).out +# Error = logs_parSub/simfit_toy_fullAngularMass_Swave_\$INT(nsamp)_\$INT(bin).err