#include #include "RooIpatia2.h" #include "RooRelBreitWigner.h" #include "RooKeysPdf.h" #include "RooChebychev.h" #include "RooWorkspace.h" using namespace RooFit ; using namespace RooStats ; int scelta = 17; void DsMuMu_TotalFit__FD() { Double_t muon_mass = 105.658; gStyle->SetOptStat(0); gStyle->SetOptTitle(0); gROOT->SetBatch(kTRUE); TFile* file = new TFile("/afs/cern.ch/work/f/fdeberna/Ds1_2460_Analysis/LAVORO/Run2__Ds1_selection/TMVA_results/ClassApp_Dsmumu/Final_BDT__Dsmumu__Run2.root"); TTree* tree = (TTree*)file->Get("Ds12DspmupmumTree"); TFile* fileWS = new TFile("13_BDT_BDTG_WS_mupm__Run2.root"); TTree* treeWS = (TTree*)fileWS->Get("Ds12DspmumuWSTree"); TH1F* HDsmumu = new TH1F("HDsmumu", "HDsmumu", 250, 2200., 2700.); tree->Draw("Ds1_M_DTF_PV>>HDsmumu"); RooRealVar Ds1_M_DTF_Ds_PV("Ds1_M_DTF_Ds_PV", "m(D_{s}^{+} #mu^{-} #mu^{+}) [MeV/c^{2}]", 2200., 2700.); //============================================================================== //============= Building PDF for Ds1(2536) ===================================== //============================================================================== RooRealVar sigmean_2536("m_{D_{s1}(2536)}", "D_{s1}(2536) mass", 2535.10, 2526., 2546.); RooRealVar sigwidth_2536("#Gamma_{D_{s1}(2536)}", "D_{s1}(2536) width2", 0.92, 1., 14.); sigmean_2536.setConstant(kTRUE); sigwidth_2536.setConstant(kTRUE); RooRealVar mDau1("mDau1", "mDau1", 1968.28); mDau1.setConstant(kTRUE); RooRealVar mDau2("mDau2", "mDau2", 2 * muon_mass); mDau2.setConstant(kTRUE); RooRealVar spin("spin", "spin", 0); spin.setConstant(kTRUE); RooRealVar radius("radius", "radius", 0.003); radius.setConstant(kTRUE); RooRelBreitWigner RBW_2536("RBW_2536", "RBW_2536", Ds1_M_DTF_Ds_PV, sigmean_2536, sigwidth_2536, spin, radius, mDau1, mDau2); TFile* file_res = new TFile("ResolutionS/Param_Fit_Mass_Resolutions.root"); RooFitResult* fitresult = (RooFitResult*)file_res->Get("ParFitResolutions"); RooRealVar* mcI1_L = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_L"); RooRealVar* mcI1_beta = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_beta"); RooRealVar* mcI1_zeta = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_zeta"); RooRealVar* mcI1_a1 = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_a1"); RooRealVar* mcI1_n1 = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_n1"); RooRealVar* mcI1_a2 = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_a2"); RooRealVar* mcI1_n2 = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_n2"); RooRealVar* dM_2536 = (RooRealVar*)fitresult->floatParsFinal().find("dM_2536"); RooRealVar* dM_2460 = (RooRealVar*)fitresult->floatParsFinal().find("dM_2460"); RooRealVar* mcI1_res_2536 = (RooRealVar*)fitresult->floatParsFinal().find("mcI1_res_2536"); RooRealVar* resratio = (RooRealVar*)fitresult->floatParsFinal().find("resratio"); RooRealVar mcI1_Lb("mcI1_Lb", "mcI1_Lb", -10.);//variabile? RooRealVar mcI1_betab("mcI1_betab", "mcI1_betab", 0.); RooRealVar mcI1_zetab("mcI1_zetab", "mcI1_zetab", 0.); RooRealVar mcI1_a1b("mcI1_a1b", "mcI1_a1b", mcI1_a1->getVal()); RooRealVar mcI1_n1b("mcI1_n1b", "mcI1_n1b", mcI1_n1->getVal()); RooRealVar mcI1_a2b("mcI1_a2b", "mcI1_a2b", mcI1_a2->getVal()); RooRealVar mcI1_n2b("mcI1_n2b", "mcI1_n2b", mcI1_n2->getVal()); RooRealVar dM_2536b("dM_2536b", "dM_2536b", 0.); RooRealVar mcI1_res_2536b("mcI1_res_2536b", "mcI1_res_2536b", mcI1_res_2536->getVal()); RooIpatia2 Ip_2536("Ip_2536", "Ip_2536", Ds1_M_DTF_Ds_PV, mcI1_Lb, mcI1_zetab, mcI1_betab, mcI1_res_2536b, dM_2536b, mcI1_a1b, mcI1_n1b, mcI1_a2b, mcI1_n2b); RooRealVar res2536("#sigma_{D_{s1}(2536)}", "D_{s1}(2536) width2", mcI1_res_2536->getVal(), 0.1, 5.); res2536.setConstant(kTRUE); Ds1_M_DTF_Ds_PV.setBins(10000, "cache"); RooFFTConvPdf PDF_Ds1_2536("PDF_Ds1_2536", "PDF_Ds1_2536", Ds1_M_DTF_Ds_PV, RBW_2536, Ip_2536); //============================================================================== //============= Building PDF for Ds1(2460) ===================================== //============================================================================== RooRealVar resratiot("resratiot", "resratiot", resratio->getVal());//variabile? RooProduct mcI1_res_2460b("mcI1_res_2460b", "mcI1_res_2460b", RooArgSet(resratiot, mcI1_res_2536b)); RooRealVar dM_2460b("dM_2460b", "dM_2460b", 0.); RooIpatia2 Ip_2460("Ip_2460", "Ip_2460", Ds1_M_DTF_Ds_PV, mcI1_Lb, mcI1_zetab, mcI1_betab, mcI1_res_2460b, dM_2460b, mcI1_a1b, mcI1_n1b, mcI1_a2b, mcI1_n2b); RooRealVar dM("m_{D_{s1}(2460)}-m_{D_{s1}(2536)}", "dM", -80., -100., -50.); RooAddition sigmean("sigmean", "dM+sigmean_2536", RooArgSet(dM, sigmean_2536)); RooRealVar sigwidth("#Gamma_{D_{s1}(2460)}", "D_{s1}(2460) width2", 0.2, 0.00000001, 10.); //If Min = 0 --> LikelihoodScan misbehaves (da Plot_..._Trigger_nsig.C) RooRelBreitWigner RBW_2460("RBW_2460", "RBW_2460", Ds1_M_DTF_Ds_PV, sigmean, sigwidth, spin, radius, mDau1, mDau2); RooFFTConvPdf PDF_Ds1_2460("PDF_Ds1_2460", "PDF_Ds1_2460", Ds1_M_DTF_Ds_PV, RBW_2460, Ip_2460); //============================================================================== //============= Building PDF for Ds1(2460)/Ds1(2536)_shifted =================== //============================================================================== TFile file_R4("/afs/cern.ch/work/f/fdeberna/Ds1_2460_Analysis/LAVORO/Run2__Ds1_selection/TMVA_results/ClassApp_2460_2536_MC/Final__2460_Dsstmumu__2016_MC.root"); TTree* tree_R4 = (TTree*)file_R4.Get("Ds12DspmupmumTree"); RooDataSet data_R4("data_R4", "dataset with x", RooArgSet(Ds1_M_DTF_Ds_PV), Import(*tree_R4)); file_R4.Close(); TFile file_R44("/afs/cern.ch/work/f/fdeberna/Ds1_2460_Analysis/LAVORO/Run2__Ds1_selection/TMVA_results/ClassApp_2460_2536_MC/Final__2536_Dsstmumu__2016_MC.root"); TTree* tree_R44 = (TTree*)file_R44.Get("Ds12DspmupmumTree"); RooDataSet data_R44("data_R44", "dataset with x", RooArgSet(Ds1_M_DTF_Ds_PV), Import(*tree_R44)); file_R44.Close(); RooKeysPdf R4_feed("R4_feed", "R4_feed", Ds1_M_DTF_Ds_PV, data_R4, RooKeysPdf::NoMirror); RooKeysPdf R44_feed("R44_feed", "R44_feed", Ds1_M_DTF_Ds_PV, data_R44, RooKeysPdf::NoMirror); //============================================================================== //============= Building Polynomial background PDF ============================= //============================================================================== RooRealVar pzero("pzero", "p0", 1.11, 0., 20.); RooRealVar pone("pone", "p1", 0.11, 0., 25.); RooChebychev poly("poly", "polinomial", Ds1_M_DTF_Ds_PV, RooArgList(pzero, pone)); ////============================================================================ RooRealVar BR_2460("BR_2460", "BR_2460", 0.1, 0., 1.); RooRealVar BR_2536("BR_2536", "BR_2536", 0.1, 0., 1.); RooRealVar nsig("nsig_{D_{s1}(2460) #rightarrow D_{s}#mu#mu}", "#signal events", 100, 0., 15000.); RooRealVar nsig2("nsig_{D_{s1}(2536) #rightarrow D_{s}#mu#mu}", "#signal events", 100, 0., 15000.); RooProduct nsig3_BR("nsig_{D_{s1}(2460) #rightarrow D_{s}*#mu#mu}", "#signal events", RooArgSet(BR_2460, nsig)); RooProduct nsig4_BR("nsig_{D_{s1}(2536) #rightarrow D_{s}*#mu#mu}", "#signal events", RooArgSet(BR_2536, nsig2)); RooRealVar nbkg("nbkg", "#background events", 2000., 0., 15000.); RooDataSet* data_DsMuMu_BR = new RooDataSet("data_DsMuMu_BR", "dataset with x", RooArgSet(Ds1_M_DTF_Ds_PV), Import(*tree)); //=================================================================================================================================================================================================================================================== //=================================================================================================================================================================================================================================================== RooRealVar nsig3("nsig_{D_{s1}(2460) #rightarrow D_{s}*#mu#mu}", "#signal events", 100, 0., 15000.); RooRealVar nsig4("nsig_{D_{s1}(2536) #rightarrow D_{s}*#mu#mu}", "#signal events", 10, 0., 15000.); RooAddPdf sum("sum", "g+g+p", RooArgList(PDF_Ds1_2460, PDF_Ds1_2536, R4_feed, R44_feed, poly), RooArgList(nsig, nsig2, nsig3, nsig4, nbkg)); // Nominal Fit RooRealVar Ds1_M23_DTF_Ds_PV("Ds1_M23_DTF_Ds_PV", "m(#mu^{-} #mu^{+}) [MeV/c^{2}]", 200., 1200.); RooRealVar Ds1_M12_DTF_Ds_PV("Ds1_M12_DTF_Ds_PV", "m(D_{s}^{+} #mu^{+}) [MeV/c^{2}]", 0., 10000.); RooDataSet* data_DsMuMu_RS2 = new RooDataSet("data_DsMuMu_RS2", "dataset with x", RooArgSet(Ds1_M_DTF_Ds_PV, Ds1_M23_DTF_Ds_PV, Ds1_M12_DTF_Ds_PV), Import(*tree)); printf("\n\nEffettuiamo il fitTo:\n\n"); RooFitResult* fitres_Dsmumu = sum.fitTo(*data_DsMuMu_RS2, Extended(), Save(kTRUE)); printf("\n\n\n\nRooWorkspace: \n\n"); RooWorkspace workspace("Dsmumu_workspace"); workspace.import(sum); workspace.import(*data_DsMuMu_RS2, Rename("Dsmumu_data")); workspace.writeToFile("WorkspaceDsMuMu.root"); workspace.Print(); }