How to get the maximum of a RooFFTConvPdf?

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