Building a pdf using RooGenericPdf and extract parameters from fitting an Asimov data

[Abdualazem]

Hello root experts,

I’m trying to build a pdf using RooGenericPdf, and then generate data using the new pdf. Hence, I would like to get the following 10 parameters: a1, a2, a3, b1, b2, b3, b4, c1, c2, c3, and c4. But I had no luck with it. Here’s what I’m trying to do:

#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooRealVar.h"
#include "RooArgSet.h"
#include "RooConstVar.h"
#include "RooProduct.h"


using namespace RooFit ;


void read(){


RooRealVar* m4l = new RooRealVar("m4l","m4l",200, 2000);

RooConstVar* m0 = new RooConstVar("m0","m0", 265);

RooConstVar a1("a1","a1", 100);
RooConstVar a2("a2","a2", 100);
RooConstVar a3("a3","a3", 100);

RooConstVar b1("b1","b1", 0.5);
RooConstVar b2("b2","b2", 0.5);
RooConstVar b3("b3","b3", 0.5);

RooConstVar c1("c1","c1", 0.0005);
RooConstVar c2("c2","c2",-0.0003);
RooConstVar c3("c3","c3", 0.00001);
RooConstVar c4("c4","c4", 0.0005);



RooFormulaVar* f1m4l = new RooFormulaVar("f1m4l","(0.5+0.5*TMath::Erf((@0-@1)/@2))*(1./(1.0+TMath::Exp((@0-@1)/@3)))", RooArgList(*m4l,a1,a2,a3));
RooFormulaVar* f1m0  = new RooFormulaVar("f1m0", "(0.5+0.5*TMath::Erf((@0-@1)/@2))*(1./(1.0+TMath::Exp((@0-@1)/@3)))", RooArgList(*m0,a1,a2,a3));

RooFormulaVar* f2m4l = new RooFormulaVar("f2m4l","TMath::Exp(@1+@2*@0)", RooArgList(*m4l,b1,b2));
RooFormulaVar* f2m0  = new RooFormulaVar("f2m0", "TMath::Exp(@1+@2*@0)", RooArgList(*m0 ,b1,b2));

RooFormulaVar* f3m4l = new RooFormulaVar("f2m4l","TMath::Exp(@1+@2*@0+@3*TMath::Power(@0,2.0)+@4*TMath::Power(@0,2))", RooArgList(*m4l,c1,c2,c3,c4));
RooFormulaVar* f3m0  = new RooFormulaVar("f2m0", "TMath::Exp(@1+@2*@0+@3*TMath::Power(@0,2.0)+@4*TMath::Power(@0,2))", RooArgList(*m0 ,c1,c2,c3,c4));

RooFormulaVar* Cnorm = new RooFormulaVar("Cnorm","@2/(@0+@1)", RooArgList(*f1m0,*f2m0,*f3m0));

RooFormulaVar* StepFuncLow  = new RooFormulaVar("StepFuncLow" ,"1.0*(@0<@1)", RooArgList(*m4l, *m0));
RooFormulaVar* StepFuncHigh = new RooFormulaVar("StepFuncHigh","1.0*(@0>=@1)",RooArgList(*m4l, *m0));

RooFormulaVar* fm4l = new RooFormulaVar("fm4l","(@0+@1)*@2*@3 + @4*@5",RooArgList(*f1m4l, *f2m4l, *StepFuncLow, *Cnorm, *f3m4l, *StepFuncHigh));

RooGenericPdf* func = new RooGenericPdf("func","@0",RooArgList(*fm4l));

f1m4l->Print();
f1m0->Print();
f2m4l->Print();
f2m0->Print();
f3m4l->Print();
f3m0->Print();
func->Print();

RooDataSet *data = func->generate(*m4l, 20000);

RooPlot* frame = m4l->frame();
func->fitTo(*data);

data->plotOn(frame);
func->plotOn(frame);
}

Any help will be very much appriciated.

Cheers,
Abdualazem.

[couet]

I think @moneta can help you.

[Abdualazem]

Hi Couet,

Thanks a lot. I modified the code as follows:

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"
#include "RooConstVar.h"

using namespace RooFit ;


void read(){
RooRealVar* m4l = new RooRealVar("m4l","m4l",140, 1800);

RooRealVar* m0 = new RooRealVar("m0","m0", 265);


RooRealVar a1("a1","a1", 4.46139e+01,0.0,4.46139e+01);
RooRealVar a2("a2","a2", 9.00000e+01,0.0,9.00000e+01);
RooRealVar a3("a3","a3", 9.00000e+01,0.0,9.00000e+01);
RooRealVar b1("b1","b1", 2.00000e-01, 0.,2.00000e-01);
RooRealVar b2("b2","b2", 1.28793e-14,0.0,1.28793e-14);
RooRealVar c1("c1","c1", 9.35895e-02,0.,9.35895e-02);
RooRealVar c2("c2","c2", 9.64895e-05,0.,9.64895e-05);
RooRealVar c3("c3","c3", 9.99999e-05,0.,9.99999e-05);
RooRealVar c4("c4","c4", 9.53674e-13,0.0,9.53674e-13);

RooFormulaVar* f1m4l = new RooFormulaVar("f1m4l","TMath::Exp(@1+@2*@0)", RooArgList(*m4l,b1,b2));
RooFormulaVar* f1m0  = new RooFormulaVar("f1m0", "TMath::Exp(@1+@2*@0)", RooArgList(*m0 ,b1,b2));

RooFormulaVar* f2m4l = new RooFormulaVar("f2m4l","(0.5+0.5*TMath::Erf((@0-@1)/@2))*(1./(1.0+TMath::Exp((@0+@1)/@3)))", RooArgList(*m4l,a1,a2,a3));
RooFormulaVar* f2m0  = new RooFormulaVar("f2m0", "(0.5+0.5*TMath::Erf((@0-@1)/@2))*(1./(1.0+TMath::Exp((@0+@1)/@3)))", RooArgList(*m0,a1,a2,a3));


RooFormulaVar* f3m4l = new RooFormulaVar("f3m4l","TMath::Exp(@1+@2*@0+@3*TMath::Power(@0,2.0)+@4*TMath::Power(@0,3))",RooArgList(*m4l,c1,c2,c3,c4));
RooFormulaVar* f3m0  = new RooFormulaVar("f3m0", "TMath::Exp(@1+@2*@0+@3*TMath::Power(@0,2.0)+@4*TMath::Power(@0,3))",RooArgList(*m0 ,c1,c2,c3,c4));

RooFormulaVar* Cnorm = new RooFormulaVar("Cnorm","@2/(@0+@1)", RooArgList(*f1m0,*f2m0,*f3m0));

RooFormulaVar* StepFuncLow  = new RooFormulaVar("StepFuncLow" ,"1.0*(@0<@1)", RooArgList(*m4l, *m0));
RooFormulaVar* StepFuncHigh = new RooFormulaVar("StepFuncHigh","1.0*(@0>@1)",RooArgList(*m4l, *m0));

RooFormulaVar* fm4l = new RooFormulaVar("fm4l","(@0+@1)*@5*@3 + @4*@2",RooArgList(*f1m4l, *f2m4l, *StepFuncLow, *Cnorm, *f3m4l, *StepFuncHigh));
RooGenericPdf* func = new RooGenericPdf("func","@0",RooArgList(*fm4l));

f1m4l->Print();
f1m0->Print();
f2m4l->Print();
f2m0->Print();
f3m4l->Print();
f3m0->Print();
func->Print();

RooDataSet *Asimvodata = func->generate(*m4l, 20000);
RooPlot* frame = m4l->frame();
func->fitTo(*Asimvodata);
Asimvodata->plotOn(frame);
func->plotOn(frame);

func->paramOn(frame,Layout(0.55)) ;
Asimvodata->statOn(frame,Layout(0.55,0.99,0.8)) ;

TText* txt = new TText(2,100,"") ;
//  txt->SetTextSize(0.04) ;
//  txt->SetTextColor(kRed) ;
  frame->addObject(txt) ;

frame->Draw() ;

}

That gives me the following plot:

Which solves the first problem I described in my first post. However, I still have a problem with the parameters as it’s not giving me what I want. I need to get a shape like the following plot:

In the code I posted above, there are three functions which are labelled f1, f2 and f3 on the last plot. I’m not sure what’s happening. It will be great If someone can help! Thanks.

Cheers,
Abdualazem.

[StephanH]

Hi Abdualazem,

to me it’s not yet clear what you want.

What exactly do you mean by “I need to get sharp”?

[Abdualazem]

Hi Stephan,

I need to get a shape similar to the second plot. But as you can see from the first plot that there’s a problem with the initial parameters. Thanks.

Cheers,
Abdualazem.

[StephanH]

Ok, so you mean that the formulas that you typed down don’t represent the distribution you want.
Unfortunately, we don’t have the time to work on such things in the root forum. It’s meant for technical questions on how to use root.

Do you have a supervisor that you could ask for ideas?
The only thing I can say at the moment is that finding starting values for parameters can be a bit tricky. Also keep in mind that the allowed range for the parameters needs to be large enough. If they hit the range limits, they will not be changed by the fitter.

The distribution you want to model looks roughly like an exponential distribution with some kind of turnon, but it would be better to think about the physics process that generates such a distribution. This should dictate the fit model.

[Abdualazem]

@StephanH

Do you have a supervisor that you could ask for ideas?

Yes, I do. I think I need to ask him.
Thanks a lot.

Cheers,
Abdualazem.