#include "RooRealVar.h" #include "RooDataSet.h" #include "RooGenericPdf.h" #include "RooArgList.h" #include "RooFitResult.h" #include "RooPlot.h" #include "RooFormulaVar.h" #include #include #include "TMinuit.h" #include #include "TFile.h" #include "TTree.h" #include "TLorentzVector.h" #include "TFitResult.h" #include "TCanvas.h" #include "TAxis.h" #include "RooClassFactory.h" #include "TROOT.h" #include "TEfficiency.h" #include #include //#include using namespace RooFit; using namespace std; void Rootfit_Corrected() { //variables globally : RooRealVar s("s", "s", 129600); RooRealVar MP("MP", "Pi constant", TMath::Pi()); RooRealVar alpha("alpha", "alpha constant", 0.0072973525643); RooRealVar mtop("mtop", "mtop", 172); RooRealVar beta_val("beta_val", "beta_val constant", 0.318433); RooRealVar MZ("MZ", "MZ", 91); RooRealVar d_p("d_p", "d_p", 0.633434); RooRealVar C("C", "C constant", 1.08); RooRealVar Bz("Bz", "Bz constant", 0.324675); RooRealVar Az("Az", "Az", 0.229703); RooRealVar By("By", "By constant", 0); RooRealVar Ay("Ay", "aY constant", -0.666667); RooRealVar B("B", "B constant", 0.516952); RooRealVar ve("ve", "ve constant", -0.08); RooRealVar CosthetaW("CosthetaW", "CosthetaW constant", 0.77); RooRealVar SinthetaW("SinthetaW", "SinthetaW constant", 0.23); RooRealVar r("r", "r constant", 0.214107); RooRealVar W("W", "W constant", 0.882105); RooRealVar Dgamma("Dgamma", "Dgamma constant", 0.05); RooRealVar x("x", "x", 0.2, 1); RooRealVar costheta("costheta", "costheta", -1, 1); RooRealVar costheta_val("costheta_val","costheta_val",-1, 1); RooRealVar x_val("x_val","x_val",0.1,1); RooRealVar E_f("E_f","E_f",0,260); // parameter to minimize : RooRealVar EVA("EVA","EVA",-0.11,-0.1,1); RooRealVar EV("EV","EV",-0.2,-0.5,1); RooRealVar EA("EA","EA",0.02,0.015,1); RooRealVar DV("DV","DV",-0.5,0.5,1); RooRealVar DVA("DVA","DVA",0.4,0.5,1); RooRealVar DA("DA","DA",0.4,0.1,1); // Omega PLus ------------------------- RooFormulaVar omega_plus( "omega_plus", "((r < B) ? ((x >= B*r && x < r) ? (1 - r) : " " (x >= r && x < B) ? (1 - x) : " " (x >= B) ? (1 - x) : 0) : 0)", RooArgList(x,r, B) ); // Omega Minus ------------------------ RooFormulaVar omega_minus( "omega_minus", "((r < B) ? ((x >= B*r && x < r) ? (1 - (x/B)) : " " (x >= r && x < B) ? (1 - (x/B)) : " " (x >= B) ? 0 : 0) : 0)", RooArgList(x,r, B) ); // F ----------------------------------------------------------------- RooFormulaVar F( "F", "((1 + beta_val)/beta_val) * (3/(2*W)) * (pow( omega_plus ,2) - pow( omega_minus ,2))", RooArgList(W, beta_val,omega_minus,omega_plus) ); // H1 ----------------------------------------------------------------- RooFormulaVar H1( "H1", "((1 - pow(beta_val, 2)) / x) * (1 / (2 * W * beta_val)) * (pow(omega_plus, 2) * (3 - 2 * omega_plus) - (pow(omega_minus, 2) * (3 - 2 * omega_minus)))", RooArgList( x,W, beta_val, omega_minus, omega_plus) ); // H2 ----------------------------------------------------------------- RooFormulaVar H2( "H2", "(1/(4*beta_val*W)) * (1 + beta_val) * pow((1 - beta_val),2) * (1/pow(x,2)) * " "(pow(omega_plus,2)*(6 - 8*omega_plus + 3*pow(omega_plus,2)) - " "pow(omega_minus,2)*(6 - 8*omega_minus + 3*pow(omega_minus,2)))", RooArgList(x,W, beta_val,omega_minus,omega_plus) ); // G ----------------------------------------------------------------- RooFormulaVar G( "G", "F + (3*pow((1 + beta_val),2)/(beta_val*W)) * x * (" "(omega_plus + std::log(std::abs(1 - omega_plus)) - " "(omega_minus + std::log(std::abs(1 - omega_minus)))))", RooArgList(x,F, W, beta_val,omega_minus,omega_plus) ); RooArgList params(EVA, EV, EA, DV, DA, DVA); RooAbsPdf *genpdf = RooClassFactory::makePdfInstance("GenPdf", " (0.5 * (3 - pow(beta_val, 2)) * (DV) - 0.5 * (1 - (3 * pow(beta_val, 2))) * (DA) - " "((1 - pow(beta_val, 2)) * (DVA)) * F + 2 * ((DVA) * G) + " "0.5 * ((DV) + (DA) + 2 * (DVA))) * (2 * H1 - H2) + " "(2 * (2 * (EVA) + (1 - pow(beta_val, 2)) * (EA)) * F + " "2 * ((EV) + (EA)) * G - 2 * (2 * (EVA) + (EV) + (EA)) * H1) * costheta + " "(0.5 * ((3 - pow(beta_val, 2)) * ((DV) + (DA)) + 6 * ((1 - pow(beta_val, 2)) * (DVA))) * F + " "2 * (DVA) * G + (-3. / 2.) * (DV + DA+ 2 * DVA) * (2 * H1 - H2)) * pow(costheta, 2)", RooArgList(x,costheta, EVA, EV, EA, DV, DA, DVA, F, G, H1, H2, beta_val)); RooDataSet* data = genpdf->generate((x,costheta),500000); RooDataHist binnedModel("binnedModel", "data", RooArgSet(x, costheta), (*data)); RooAbsReal* nll = genpdf->createNLL(binnedModel, RooFit::Extended(true), RooFit::NumCPU(4)); RooMinimizer minim(*nll); minim.setStrategy(2); minim.setEps(1); minim.minimize("Minuit2", "scan"); minim.minimize("Minuit2", "migrad"); minim.minimize("Minuit2", "hesse"); RooPlot* frameEA = EA.frame( RooFit::Range(0.2, 1.0), RooFit::Title("NLL vs EA")); frameDA->SetXTitle("EA"); frameDA->SetYTitle("NLL"); // best-fit values double bestEA = EA.getVal(); double nll_min = nll->getVal(); RooAbsReal* pll_EA = nll->createProfile(EA); pll_EA->plotOn(frameEA, RooFit::LineColor(kBlue), RooFit::Precision(1e-4), RooFit::DrawOption("L") ); TCanvas* c = new TCanvas("c", "NLL vs EA", 800, 600); c->SetLeftMargin(0.15); c->SetBottomMargin(0.15); frameEA->Draw(); frameEA->GetYaxis()->SetTitleOffset(1.4); gPad->SetGridx(); gPad->SetGridy(); }