Hi,
I have a 1D histogram to which I fit a RooAddPdf (which is 2 gaussians+ a straight line).
I then want to know, for a certain value of the variable, the height/value of the RooAddPdf.
There is something in the normal TF1, called Eval(x) and that give you value of the TF1 at x.
But for a RooAddPdf, I cannot find any function that does it.
Can someone point me to the right method?
Regards
Fatima Soomro
Hello,
You can use RooAbsPdf::getVal()
see: root.cern.ch/root/html/RooAbsPdf … Pdf:getVal
Regards,
– Gregory
Hi,
thanks for that.
Now I have another question:
I am fitting a gaussian+chebychev + a mass distribution, and when i get the integral of the pdfs (using getNorm)
norm: signal 6.12803 bgk 1 total 1
how is it that the total pdf integral < signal integral+bkg integral
am I missing something here?
And, I have been trying to follow the documentation on accessing a ttree and filling a RooDataSet.
Unfortunately, without any luck… On the webpage:
RooDataSet(const char* name, const char* title, TTree* t, const RooArgSet& vars, const RooFormulaVar& cutVar, const char* wgtVarName = 0)
can someone elaborate what is the RooArgSet vars is here?
And also the RooFormulaVar cutVar? can I use a string here?
I have a TTree object called t, and a cut string called cut. I want to plot variable “m” while making some cuts on some variables in t so I have tried:
RooDataSet data(“data”, “data”, t,RooArgSet(m),cut);
Regards
Fatima Soomro
Hi Fatima,
The normalization integral that is returned by getNorm() is not the integral over the pdf, it is the term by wich getVal() needs to be divided to obtain a normalized pdf. Every normalized pdf integrated over its
nominal observable range will yield exactly one (by definition). Alternatively you can integrate a pdf
over a subrange of its observables e.g.
x.setRange(“A”,5,10) ;
RooAbsReal* intPdf = pdf.createIntegral(x,x,“A”) ;
which will return the fraction of the pdf contained in range “A”, which will be a number between 0 and 1 provided range “A” is contained in the nominal range of the observables.
Concerning your dataset question: If you want to create a dataset from a tree you first have to
define the observables as RooRealVar objects. E.g. for a case with two observables you an do
RooRealVar x(“x”,“x”,-10,10) ;
RooRealVar y(“y”,“y”,-10,10) ;
RooDataSet d(“d”,“d”,RooArgSet(x,y),Import(*t)) ;
which will import the data contained in the tree branches “x” and “y” contained in the given tree t. Please note that at the moment you can only make cuts on variables you import, so if you add a Cut(“y<0”) it
would work, but if you want to cut on another (hypothetical) observable “z”, you’d have to import that one too. [ This is a technical limitation that will be removed some time this year ]
Wouter
I would like to resurrect the original question. Based on the documentation, I’m unclear on how to use getVal() to obtain the height/value of the pdf at a specific point. Thus far my attempts have yielded very inaccurate results. As an example:
[code]//retrieve normalized value:
RooRealVar myTemp(“myTemp”,“myTemp”,490,510); //interested in height/value at 500
RooArgSet nset(myTemp);
double normExpVal = myPDF.getVal(&nset);
//plot constant value (check that it intercepts myPDF at 500):
RooRealVar con(“con”,“con”,normExpVal,normExpVal,normExpVal);
con.plotOn(xframe,LineColor(kViolet));[/code]
Please see the attached image for the results. Ideally, the violet line should intersect the red curve at 500 on the x-axis, but this is clearly not the case. In fact, in two of the cases the value is so ridiculous that the remaining curves (the code for which is not shown here) aren’t even visible, while in the first case it’s too small to be shown itself. The precise values, for those who are curious, are these:
normExpVal 1 = 0.00588492
normExpVal 2 = 0.0164406
normExpVal 3 = 0.0121013
normExpVal 4 = 1.68935e+39
normExpVal 5 = 1.8627e+39
normExpVal 6 = 0.0180288
Could someone shed some light on the drastic range of these results, or ideally, how to use getVal() to retrieve the height at a specific point?
Many thanks,
Robert
Hi, thanks for resurrecting it
I have also been unable to understand what is being returned from my fits.
So I construct my pdf from a gaussain and a chebychev
RooAddPdf model(“model”,“model”, RooArgList(gauss, bkg), fsig);
in a variable m, which is a rooRealVar
Then I fit my data points with this pdf, and to find the value of the pdf at a certain value, call it phim, of my variable m, I do
m.setVal(phim);
gauss.evaluate();
model.evaluate();
fsig.getVal();
bkg.evaluate();
For the attached plot, I get the following values for
phim= 1019.11 (which is very signal like)
gauss: 0.961827 ,fsig: 0.187641, bkg: 1.26698, model: 1.20972
so weight of this candidate, with mass 1019 is
(gauss*fsig)/model= 0.149191
which does not look right to me, I think it should have a higher weight
Regards
Fatima Soomro
Hi Fatima,
The “model” I’m using is actually a RooGenericPdf, and my attempts to duplicate your approach were met with errors when I tried to call evaluate() on such a model (though I do not know why). I did however try setting the variable via setVal(phim) prior to calling getVal(), but the results were similarly inaccurate.
Does anyone know the difference between evaluate() and getVal(), and which is more appropriate to this application?
Thanks,
Robert
Hi,
I’m not if this answer all questions, but I can clarify a few issues here.
A RooFit pdf does not have an intrinsic definition of how it is normalized as a pdf
(there are good reasons for this on which I will not elaborate now). So if you want
to have a normalized probability density from a RooFit pdf you always need to pass
it a set of variable that should be interpreted as observables, e.g. for a Gaussian
// Assume a RooGaussian g(x,m,s)
RooArgSet obs(x) ;
Double_t val = gauss.getVal(&obs) ;
returns the normalized probability density. Merely calling
Double_t val = gauss.getVal()
will return the value of RooGaussian::evaluate() which is not guaranteed to be normalized properly.
Note that you can also call
RooArgSet obs2(x,s) ;
Double_t val = gauss.getVal(&obs2) ;
which would interpret gauss as a 2D pdf with observables x and s. Generally F.getVal(obs) returns
F(x,obs) = f(x) / Int d(obs) f(x)
where f(x) is the value returned by evaluate(). The normalization integral is an analytical expression for many pdfs (like RooGaussian). In case no analytical expression is available (as e.g. for RooGenericpdf) a numeric integral is used.
Wouter
Hi Folks,
This’s still about the original question. Here’s my set-up:
// missing mass2 for B -> rho l (nu)
RooRealVar mm2(“mm2”,“mm2”,0,-0.5,2);
// some signal PDF
RooHistPdf sighistpdf(“sighistpdf”,“sighistpdf”,mm2,roohistpdf,2) ;
// some background function ansatz
RooDoubGausPdf bkg(“bkg”,“background p.d.f.”,mm2,m0,s0,c1,m1,s1) ;
// construct the model where nsig and nbkgd will be returned by the fit
RooAddPdf model(“model”,“model”,RooArgList(sighistpdf,bkg),RooArgList(nsig,nbkgd));
// do an extended ML fit …
// I can plot the signal component using
model.plotOn(xframe,Components(sighistpdf),LineColor(2),LineStyle(kDashed));
Now, what’s the simplest way I can get the prediction of the signal fraction from the model, at a particular value of mm2? For example, I can do
mm2.setVal(0);
then what?
Thanks!
Furthernore…
I tested with
TH1F *h2 = new TH1F(“h2”,"",50,-0.5,2.0);
for(float x = -0.5; x < 2.0; x += dx){
mm2.setVal(x);
h2->Fill(x,sighistpdf.getVal());
}
h2 certainly gives the correct shape. The problem is with the normalization, which depends on the binning, range, etc. I mean, when I’m plotting the components, there should be a simple way to ask what is the value of the ratio: sig(mm2_i) / ( sig(mm2_i) + bkgd(mm2_i) ) at a particular point mm2_i.