Include all changes since unblinding

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)

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; ParIndx<nPars; ++ParIndx){
+      RooRealVar* Par = (RooRealVar*)fitResult->floatParsFinal().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; ParIndx<nPars; ++ParIndx){
+    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][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<string> files, std::vector<string> labels, bool datalike)
+{
+
+  double pVal [nPars];
+  double pErrH [nPars];
+  double pErrL [nPars];
+
+  for (int i=0; i<nPars; ++i) {
+    pVal [i] = 0.5 + 1*i;
+    pErrH[i] = 0.5 ;
+    pErrL[i] = 0.5 ;
+  }
+
+  double genRes [nPars];
+  double genErrH [nPars];
+  double genErrL [nPars];
+
+  unsigned int n_files = files.size();
+  string stat = datalike ? "_dataStat":"_MCStat";
+
+  double Res    [n_files][nPars];
+  double Diff    [n_files][nPars];
+  double ErrH    [n_files][nPars];
+  double ErrL    [n_files][nPars];
+
+  double DiffErrH    [n_files][nPars];
+  double DiffErrL    [n_files][nPars];
+
+  double pValReco  [n_files][nPars];
+  double pErrHReco [nBins];
+  double pErrLReco [nBins];
+
+  double deltaQ = 0.5;
+  for (int i=0; i<nPars; ++i) {
+    deltaQ = (1.) / n_files;
+    for (unsigned int ifile = 0; ifile < n_files; ifile++) {
+      pValReco[ifile][i] = deltaQ/2 + deltaQ*ifile + i;
+    }
+      pErrHReco[i] = 0.5 * deltaQ;
+      pErrLReco[i] = 0.5 * deltaQ;
+  }
+
+  std::vector<TGraphAsymmErrors*> 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"<<endl;
+    return;
+  }
+  RooWorkspace* wspRes = (RooWorkspace*)finGen->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<string> 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);
+
+}