void ProblemMacro () { gROOT->SetBatch(kTRUE); gROOT->SetStyle("Plain"); gSystem->Load("libRooFit.so"); using namespace RooFit; const double m_pion = 139.57018; // MeV/c^2 const double m_kaon = 493.677; // MeV/c^2 RooNumIntConfig& conf = RooNumIntConfig::defaultConfig(); conf.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); // Observables declaration RooRealVar *mpipi = new RooRealVar("mpipi","mpipi",4900,5700,"MeV/c^{2}"); RooRealVar *pidp = new RooRealVar("pidp","pion consistency of first track",-500,500); RooRealVar *pidm = new RooRealVar("pidm","pion consistency of second track",-500,500); RooRealVar *pA = new RooRealVar("pA","momentum asymmetry",-0.9,0.9,"MeV/c^{2}"); RooArgSet* pdf_obs = new RooArgSet(*mpipi,*pA,*pidp,*pidm); //Initialize parameters RooRealVar* mB0 = new RooRealVar("mB0_mean","B0 mass",5280,"MeV/c^{2}"); RooRealVar* smB1 = new RooRealVar("smB1","B0 mass resolution",20,"MeV/c^{2}"); RooRealVar* smB2 = new RooRealVar("smB2","B0 mass resolution",120,"MeV/c^{2}"); RooRealVar* normVar = new RooRealVar("normVar","Normalization",0.80); //Initialize the PDFs parameters RooRealVar* pAp1 = new RooRealVar("pAp1","pAp1",-0.0262); RooRealVar* B02Kpi_A = new RooRealVar("B02Kpi_A","B0 -> Kpi charge asymmetry",+0.008); RooRealVar* B02pipi_B02Kpi = new RooRealVar("B02pipi_B02Kpi","Br(B0->pipi)/Br(B0->Kpi)",0.27); RooRealVar* B02KK_B02Kpi = new RooRealVar("B02KK_B02Kpi","Br(B0->KK)/Br(B0->Kpi)",0.); // take 90% UL from PDG RooRealVar *N_B02pipiSig = new RooRealVar ("N_B02pipiSig","Signal number B0 -> pi+ pi-",20000); RooRealVar* N_B02KpiSig = new RooRealVar ("N_B02KpiSig","Signal number B0 -> K+ pi-",18000); pidp->setRange("pipi",pidp->getMin(),0); pidp->setRange("piK",pidp->getMin(),0); pidm->setRange("pipi",pidm->getMin(),0); pidm->setRange("piK",0,pidm->getMax()); pidp->setRange("Kpi",0,pidp->getMax()); pidp->setRange("KK",0,pidp->getMax()); pidm->setRange("Kpi",pidm->getMin(),0); pidm->setRange("KK",0,pidm->getMax()); RooAbsPdf* pdf_pidpip = new RooGaussian("pdf_pidpip","Pid p+",*pidp,RooConst(-40),RooConst(0.1)); RooAbsPdf* pdf_pidpim = new RooGaussian("pdf_pidpim","Pid p-",*pidm,RooConst(-40),RooConst(0.1)); RooAbsPdf* pdf_pidKp = new RooGaussian("pdf_pidKp","t pid for P+",*pidp,RooConst(40),RooConst(0.1)); RooAbsPdf* pdf_pidKm = new RooGaussian("pdf_pidKm","t pid for P-",*pidm,RooConst(40),RooConst(0.1)); RooAbsPdf* pdf_MomAsym = new RooGenericPdf("pdf_MomAsym", "Momenta Asymmetry PDF", "1/(1+exp((@0+1)**@1)+exp((-@0+1)**@1))", RooArgList(*pA,*pAp1)); //Unextend Mass B->ff pdf RooAbsReal* N_B02pipi = new RooFormulaVar("N_B02pipi","B0 -> pi+pi- yield","@0",RooArgList(*N_B02pipiSig)); RooGaussian *pdf_core_B02pipi = new RooGaussian("pdf_core_B02pipi","core B02pipi",*mpipi,*mB0,*smB1); RooGaussian *pdf_tail_B02pipi = new RooGaussian("pdf_tail_B02pipi","tail B02pipi",*mpipi,*mB0,*smB2); RooAbsPdf* pdf_unm_B02pipi = new RooAddPdf("pdf_sunm_B02pipi","B0->pi+pi- signal mass", *pdf_core_B02pipi, *pdf_tail_B02pipi, *normVar); //Extended Mass PDFs RooAbsPdf* pdf_m_B02pipi = new RooExtendPdf("pdf_m_B02pipi","B0 -> pi+pi- mass sig", *pdf_unm_B02pipi, *N_B02pipi); RooAbsPdf* tpdf_cnd_B02pipi = new RooProdPdf("pdf_cnd_B02pipi","conditional B0->pi+pi- signal mass", RooArgSet(*pdf_MomAsym),Conditional(RooArgSet(*pdf_m_B02pipi),RooArgSet(*mpipi)),Conditional(RooArgSet(*pdf_pidpip),RooArgSet(*pidp)),Conditional(RooArgSet(*pdf_pidpim),RooArgSet(*pidm))); // derived parameters RooAbsReal* N_B02Kpi = new RooFormulaVar("N_B02Kpi", "B0 -> K+pi- yield" ,"0.5*(1-@0)*@1", // Changed convention respect to v5r2 version RooArgList(*B02Kpi_A,*N_B02KpiSig)); RooAbsReal &d2 = RooRealConstant::value(m_kaon*m_kaon - m_pion*m_pion); RooAbsReal *mB0Kpi = new RooFormulaVar("mB02Kpi","sqrt(@0*@0-@1*2/(1+@2))", RooArgList(*mB0,d2,*pA)); RooGaussian *pdf_core_B02Kpi = new RooGaussian("pdf_core_B02Kpi","core B02pipi",*mpipi,*mB0Kpi,*smB1); RooGaussian *pdf_tail_B02Kpi = new RooGaussian("pdf_tail_B02Kpi","tail B02pipi",*mpipi,*mB0Kpi,*smB2); RooAbsPdf *pdf_unm_B02Kpi = new RooAddPdf("pdf_sunm_B02Kpi", "B0->K+pi- signal mass", *pdf_core_B02Kpi, *pdf_tail_B02Kpi, *normVar); RooAbsPdf *pdf_m_B02Kpi = new RooExtendPdf("pdf_m_B02Kpi", "B0 -> K+pi- mass sig", *pdf_unm_B02Kpi, *N_B02Kpi ); RooAbsPdf *tpdf_cnd_B02Kpi = new RooProdPdf("pdf_cnd_B02Kpi", "conditional B0->K+pi- signal mass", RooArgSet(*pdf_MomAsym),Conditional(RooArgSet(*pdf_m_B02Kpi),RooArgSet(*mpipi)),Conditional(RooArgSet(*pdf_pidKp),RooArgSet(*pidp)),Conditional(RooArgSet(*pdf_pidpim),RooArgSet(*pidm))); RooAbsReal* N_B02piK = new RooFormulaVar("N_B02piK", "B0 -> pi+K- yield" ,"0.5*(1+@0)*@1", // Changed convention respect to v5r2 version RooArgList(*B02Kpi_A,*N_B02KpiSig)); RooAbsReal *mB0piK = new RooFormulaVar("mB02piK","sqrt(@0*@0-@1*2/(1-@2))", RooArgList(*mB0,d2,*pA)); RooGaussian *pdf_core_B02piK = new RooGaussian("pdf_core_B02piK","core B02pipi",*mpipi,*mB0piK,*smB1); RooGaussian *pdf_tail_B02piK = new RooGaussian("pdf_tail_B02piK","tail B02pipi",*mpipi,*mB0piK,*smB2); RooAbsPdf* pdf_unm_B02piK = new RooAddPdf("pdf_sunm_B02piK", "B0->pi+K- signal mass", *pdf_core_B02piK, *pdf_tail_B02piK, *normVar); RooAbsPdf* pdf_m_B02piK = new RooExtendPdf("pdf_m_B02piK", "B0 -> pi+K- mass sig", *pdf_unm_B02piK, *N_B02piK ); RooAbsPdf* tpdf_cnd_B02piK = new RooProdPdf("pdf_cnd_B02piK", "conditional B0->pi+K- signal mass", RooArgSet(*pdf_MomAsym),Conditional(RooArgSet(*pdf_m_B02piK),RooArgSet(*mpipi)),Conditional(RooArgSet(*pdf_pidpip),RooArgSet(*pidp)),Conditional(RooArgSet(*pdf_pidKm),RooArgSet(*pidm))); RooArgList FitPdfs; FitPdfs.add(*tpdf_cnd_B02pipi); FitPdfs.add(*tpdf_cnd_B02piK); FitPdfs.add(*tpdf_cnd_B02Kpi); RooAbsPdf* pdf_FitPdf = new RooAddPdf("pdf","Total Rate with background",FitPdfs); double exp_evt = pdf_FitPdf->expectedEvents(*pdf_obs); int f_nEvt = (int)exp_evt; cout<<"Generating n. events: "<generate(*pdf_obs,f_nEvt); // Variation prior mean for simulation of multiple experiment RooNLLVar* nll_fit = new RooNLLVar("nll_fit","-log(L)",*pdf_FitPdf,*data_set,Extended(true)); RooRealVar * offset = new RooRealVar("offset","offset",0); *offset = nll_fit->getVal()+10; RooFormulaVar* L= new RooFormulaVar("L","@0-@1",RooArgList(*nll_fit,*offset)); RooMinuit m1(*L); m1.setVerbose(kTRUE) ; m1.setLogFile("myLog.txt"); m1.setProfile(1); m1.setPrintLevel(3); m1.setStrategy(2); B02Kpi_A->setConstant(kFALSE); N_B02KpiSig->setConstant(kFALSE); N_B02pipiSig->setConstant(kFALSE); mB0->setConstant(kFALSE); normVar->setConstant(kFALSE); pAp1->setConstant(kFALSE); smB1->setConstant(kFALSE); smB2->setConstant(kFALSE); res = m1.save(); res->Print("v") ; m1.migrad() ; res = m1.save(); std::cout << "fit result on NLL" << std::endl ; res->Print("v") ; res->Print() ; TCanvas *c4 = new TCanvas(); c4->cd(); RooPlot *p3 = mpipi->frame(); data_set->plotOn(p3); pdf_FitPdf->plotOn(p3); pdf_FitPdf->plotOn(p3,Components("*B02pipi*"),LineStyle(kDashed), LineColor(46)); pdf_FitPdf->plotOn(p3,Components("*B02Kpi*"), LineColor(kRed)); pdf_FitPdf->plotOn(p3,Components("*B02piK*"), LineColor(kBlack)); p3->Draw(); c4->Print("mass.eps"); return; }