Dear all,
I’m new using roofit, I need to do fits of convolution
of a Landau and a Gaussian distribution, so I try the
example “rf208_convolution”, but I would like to have
the maximum peak position and not the one which is
shifted due the Landau definition.
Is there a method to get the maximum of the pdf as
getMaximum of a TF1 function? or how could I can
get the Maximum value (and error) of the final
convoluted pdf?
Thanks for help,
Andrea
Hi,
There is currently no such function in RooAbsPdf, but looking at what is in TF1 now, I think it is easy to add that. I’ll try to get that in for the next release.
Wouter
Has such a function already been added?
I found nothing in 5.27.02 
How could one proceed in that case?
Cheers, Z
Hi,
you don;t really need such method, you can easily convert a RooAbsPdf in a TF1, just do:
RooRealVar * x = ..... // your observable x
RooAbsPdf * pdf = ..... // your pdf
TF1 * f = pdf->asTF( RooArgList(*x) );
double xmax = f->GetMaximumX(); Lorenzo
Grazie Lorenzo.
I should have looked to the ROOT page documentation :
root.cern.ch/root/htmldoc/RooAbsReal.html
instead of the sourceforge one which is not uptodate
roofit.sourceforge.net/docs/clas … sReal.html
Cheers,
Z
BTW, either in interpreted or compiled mode, one has to clone the function to be able to draw it!
(seg. fault otherwise)
RooRealVar x ("x" ,"x" ,-10.,10.) ;
RooRealVar mean ("mean" ,"mean" , 1.) ;
RooRealVar sigma ("sigma","sigma", 2.) ;
RooGaussian roopdf_gaus ("roopdf_gaus","roopdf_gaus",x,mean,sigma) ;
TF1 * f_roopdf_gaus = roopdf_gaus.asTF (RooArgList(x)) ;
//f_roopdf_gaus->Draw() ; /// *** Break *** segmentation violation
TF1 * fct = (TF1*) f_roopdf_gaus->Clone("fct") ;
fct->Draw() ;
Cheers, Z
You need to maintain the underlined pdf object alive, i.e. do
RooGaussian * roopdf_gaus = new RooGaussian("roopdf_gaus","roopdf_gaus",x,mean,sigma) ;or you can do a TF1::Clone or use TF1::DrawClone()
Lorenzo
The heap solution (RooGaussian*) does not work
(Using TF1::DrawClone() works in any cases of course)
Z
BTW, how could one get the TF1 as it is drawn on a RooPlot?
I just can’t get the same function… (different normalization)
[code]{
Double_t xmin=-10., xmax=10. ;
TH1D * histo = new TH1D (“histo”,“histo”,50,xmin,xmax) ;
histo->FillRandom(“gaus”,1000) ;
RooRealVar x (“x”,“x”,xmin, xmax) ; x.Print() ;
RooDataHist data (“data”,“data”,x,histo) ; data.Print() ;
RooRealVar mean (“mean” ,“mean” , 0., xmin, xmax) ; mean .Print() ;
RooRealVar sigma (“sigma”,“sigma”, 1., xmin, xmax) ; sigma.Print() ;
RooGaussModel roopdf_gaus (“roopdf_gaus” ,“roopdf_gaus” ,x,mean,sigma) ; roopdf_gaus.Print() ;
roopdf_gaus.fitTo(data,"") ;
TF1 * f_roopdf_gaus = roopdf_gaus.asTF(RooArgList(x)) ;
RooPlot * frame = x.frame() ;
data .plotOn(frame) ;
roopdf_gaus.plotOn(frame) ;
frame->Draw() ;
TCanvas * canvas = new TCanvas (“canvas”,“canvas”) ;
histo->Scale(1./histo->Integral()) ;
histo->Draw() ;
f_roopdf_gaus->DrawClone(“same”) ;
}
[/code]Cheers,
Z
Hi,
you have to divide the histogram also by the bin width. Do
histo->Scale(1./histo->Integral(),"width");Cheers, Lorenzo
Grazie Lorenzo 
Z
Is it possible to get back a non-normalized TF1 that can be drawn on the original data (like with the RooPlot)?
I did not manage to scale the returned TF1.
TIA,
Z
Since, you cannot easily scale the TF1, you can however create the TF1 from a scaled Pdf.
A possible way to do this is in your code is by adding these lines:
double scale_factor = histo->Integral()*histo->GetBinWidth(1);
TString formula = TString::Format("%f * roopdf_gaus",scale_factor);
// create a scaled Gaussian function = scale * Gaus
RooFormulaVar scaled_gaus("scaled_gaus",formula,RooArgList(roopdf_gaus) );
TF1 * f_roopdf_gaus = scaled_gaus.asTF(RooArgList(x) ) ;and now your TF1 will be normalized to the data, i.e you can plot directly on the histogram as it is done in the RooPlot
Lorenzo
Thanks !!!
Z
Unfortunately, either methods (histogram or roopdf scaling) does not seem to work with a Roo*ConvPdf
Double_t xmin=-5.5, xmax=20.5 ;
TH1D *histo = new TH1D("histo","histo",26,xmin,xmax);
histo->SetBinContent(1,6.764028e-07);
histo->SetBinContent(2,7.382022e-05);
histo->SetBinContent(3,0.003090511);
histo->SetBinContent(4,0.05065845);
histo->SetBinContent(5,0.337311);
histo->SetBinContent(6,0.9789554);
histo->SetBinContent(7,1.41928);
histo->SetBinContent(8,1.295449);
histo->SetBinContent(9,0.9710827);
histo->SetBinContent(10,0.6984009);
histo->SetBinContent(11,0.5004868);
histo->SetBinContent(12,0.358615);
histo->SetBinContent(13,0.2569589);
histo->SetBinContent(14,0.1841191);
histo->SetBinContent(15,0.1319271);
histo->SetBinContent(16,0.09452989);
histo->SetBinContent(17,0.06773362);
histo->SetBinContent(18,0.04853326);
histo->SetBinContent(19,0.0347756);
histo->SetBinContent(20,0.02491781);
histo->SetBinContent(21,0.01785439);
histo->SetBinContent(22,0.01279313);
histo->SetBinContent(23,0.009162816);
histo->SetBinContent(24,0.006503793);
histo->SetBinContent(25,0.004277093);
histo->SetBinContent(26,0.0021264);
/// Observables
RooRealVar * t = new RooRealVar ("t","t",xmin,xmax) ; t->Print() ;
/// Landau (t,ml,sl) ;
RooRealVar * ml = new RooRealVar ("ml","mean landau" ,0.,xmin,xmax) ; ml->Print() ;
RooRealVar * sl = new RooRealVar ("sl","sigma landau" ,1.,xmin,xmax) ; sl->Print() ;
RooLandau * roopdf_landau = new RooLandau ("roopdf_landau","roopdf_landau",*t,*ml,*sl) ; roopdf_landau->Print() ;
/// Gaussian (t,mg,sg)
RooRealVar * mg = new RooRealVar ("mg","mean gaussian" ,1.,xmin,xmax) ; mg->Print() ;
RooRealVar * sg = new RooRealVar ("sg","sigma gaussian",1.,xmin,xmax) ; sg->Print() ;
RooGaussian * roopdf_gauss = new RooGaussian ("roopdf_gauss" ,"roopdf_gauss" ,*t,*mg,*sg) ; roopdf_gauss->Print() ;
/// Convolution
t->setBins(10000,"cache") ;
RooFFTConvPdf * roopdf_landgaus_conv = new RooFFTConvPdf ("roopdf_landgaus_conv","roopdf_landgaus_conv",*t,*roopdf_landau,*roopdf_gauss) ; roopdf_landgaus_conv->Print() ;
/// RooData
RooDataHist * data = new RooDataHist ("data","data",RooArgSet(*t),histo) ; data->Print() ;
/// Fit
RooFitResult * fitresult = roopdf_landgaus_conv->fitTo (*data) ;
/// Rooplot
RooPlot * frame = t->frame() ;
data ->plotOn(frame) ;
roopdf_landgaus_conv->plotOn(frame) ;
frame->Draw() ;
/// TCanvas
TF1 * f_roopdf_landgaus_conv = roopdf_landgaus_conv->asTF(RooArgList(*t)) ;
printf ("\nf_roopdf_landgaus_conv->GetMaximum()=%g\n",f_roopdf_landgaus_conv->GetMaximum()) ;
printf ("\nhisto->GetMaximum()=%g\n",histo->GetMaximum()) ;
histo->Scale(1./histo->Integral(),"width") ;
printf ("\nhisto->GetMaximum()=%g\n",histo->GetMaximum()) ;
TCanvas * canvas = new TCanvas ("canvas","canvas") ;
canvas->Divide(1,2) ;
canvas->cd(1) ; histo->Draw("hist") ;
canvas->cd(2) ; f_roopdf_landgaus_conv->DrawClone() ;
Any clue?
Cheers,
Z
Please also mind the strange shape of the landgaus queue (from 16.5).
Using RooNumConvPdf yields a stiffer fall appearing later though (from 19.
.
Any tips to avoid this?
Hi,
in ROOT 5.27/04 (out this week) you can do the following
RooAbsPdf* myFunc ;
// Create function that represents derivative
RooAbsReal* deriv = myFunc->derivative(obs,1) ;
// Find point where derivative is zero
Double_t xmax = deriv.findRoot(obs,xmin,xmax,0) ;
Both functions already existed in earlier releases, but there was a bug
in findRoot (it returned the status code rather than the found value)
Wouter
Thanks Wouter.
But I can only get the maximum ordinate of the “not superimposable” pdf (not the “plotOn()” one),
i.e. with the previous piece of code
// roopdf derivative
RooAbsReal * roopdf_landgaus_conv_deriv = (RooAbsReal*) roopdf_landgaus_conv->derivative(*t,1) ;
// maximum abcissa
Double_t x_ymax = roopdf_landgaus_conv_deriv->findRoot(*t,xmin,xmax,0) ;
// corresponding maximum ordinate
Double_t ymax = roopdf_landgaus_conv->getVal(new RooArgSet(RooRealVar("roo_x_ymax","roo_x_ymax",x_ymax))) ;
printf ("\nx_ymax=%g / ymax=%g\n",x_ymax,ymax) ; /// returns ymax=315.392 instead of something greater than histo->GetMaximum()=1.41928 (see plot)
What I’d basically like to do is to reproduce the Rooplot on a TCanvas.
How could I do that?
Cheers,
Z
Hi,
You get an unnormalized return value in your latest code snippet because the
names of the normalization observables you pass don’t match those of the
pdf.
Generally,given a pdf F(x,y) that you want to evaluate at x=x0 and y=0 you should
do
x.setVal(x0) ;
y.setVal(y0) ;
cout << F.getVal() << endl ; // returns unnormalized value of F(x,y) ;
cout << F.getVal(RooArgSet(x,y)) << endl ; // returns F(x,y)/Int F(x,y) dx dy
what that modification, it should work.
Wouter