Plotting a RooFormulaVar

Roofit and RooStats
1

Hello,

I am using roofit and roostats to perform some analyses.
I am using ROOT 5.27/01 (trunk@32073, Jan 25 2010, 21:29:56 on linux)

I would like to make a plot of a RooFormulaVar. How does one do this?

For example,
I define some variables:

  RooRealVar x("x",  "my x var", 50.0, 1.0, 100);  
  RooRealVar y("y",  "my y var", 25.0, 1.0, 100);

Then I set up my model and import my data, then perform a fit:

result = model->fitTo(*data, Extended(kTRUE),SumW2Error(kTRUE), PrintLevel(1), Range("cutrange"), Save());

Now I want to make a plot of the ratio x/y.

    RooFormulaVar ratio_func("ratio","@0/@1",RooArgList(*x,*y));
    RooRealVar* ratio = (RooRealVar*) data->addColumn(ratio_func) ;  // line 184

    RooPlot* ratio_frame = ratio->frame(Bins(10), Title("My Ratio")) ;
    ratio->plotOn(ratio_frame, LineColor(kBlue));
    ratio_frame->Draw() ;

The code above doesn’t work. Any ideas why?
Here’s what happens:

[code] *** Break *** segmentation violation
The lines below might hint at the cause of the crash.
If they do not help you then please submit a bug report at
http://root.cern.ch/bugs. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.

#10 0x0340b5cd in std::ostream::sentry::sentry () from /usr/lib/libstdc++.so.6
#11 0x0340e00a in std::ostream::operator<< () from /usr/lib/libstdc++.so.6
#12 0x0171844e in RooDataHist::printValue () from /home/palladino/root_dev/lib/libRooFitCore.so
#13 0x0807258f in getML (ws=0x9428ec8) at getlikelihood.C:184

[/code]

Note: I can NOT redefine

  RooRealVar x("x",  "my x var", 50.0, 1.0, 100); 
  RooRealVar xOvery("xOvery","my ratio", 50.0);
  RooFormulaVar y("y","my y var", "@0/@1", RooArgList(x,xOvery)) ;

because x and y are highly correlated.

Eventually I would like to get a profile likelihood, with both x and y as parameters of interest,

I can plot the 2-D 1,2,3 sigma contours of the nll fine, but I’d really like to plot x/y of the pll.

Many thanks in advance.[/code]

2

Hi,

What you do should work. Can you see if the following modification
helps

RooPlot* ratio_frame = ratio->frame(Bins(10), Title("My Ratio"), Range(0,10)) ;

i.e. specify an explicit range for the ratio plot frame

Wouter

3

Hi Wouter,

The segmentation fault happens at line 184 (the addColumn function call)

So including a plotting range would not help. Any other suggestions?

Thanks,
Anthony

4

Hi,

Sorry, I misunderstood. I will look into this and get back to you.

Wouter

5

Hi,

I just tried this myself, but I don’t manage to reproduce the issue:

// Make dummy dataset D(x,y)
RooWorkspace w(“w”,1) ;
w.factory(“PROD::gxy(Gaussian::gx(x[1,10],5,3),Gaussian::gy(y[1,10],5,3))”) ;
RooDataSet* d = w::gxy.generate(RooArgSet(w::x,w::y),10000) ;

// Make ratio function and add to data
RooFormulaVar rfunc(“r”,"@0/@1",RooArgList(w::x,w::y)) ;
d->addColumn(rfunc) ;

// Make plot of ratio distribution
w.import(d) ;
RooPlot frame = w::r.frame(0,10) ;
d->plotOn(frame) ;
frame->Draw() ;

So it might be something specific to your data. If your data file is not too large, I can have another look at this if you send it to me.

Wouter

6

Hi Wouter,

You are absolutely correct. There was a “problem” with my data. I was trying to addColumn after I had already binned my dataset. Calling addColumn to a binned dataset resulted in the segmentation fault.

I have another question, also about plotting a different ratio x/y:

I have a RooProfileLL with two parameters of interest.

RooProfileLL pll("pll", "", *nll, RooArgSet(*y, *x));

I can plot y versus x with no problems.

TH1* h2d = pll.createHistogram("h2d",*x,Binning(num_bins),YVar(*y,Binning(num_bins)),Scaling(kFALSE)) ;

What I would really like is x/y.

RooProfileLL pllratio("pllratio", "", *nll, RooArgSet(*ratio)); TH1* hratio = pllbr->createHistogram("hratio",*ratio,Binning(num_bins_ratio),Range(0,10),Scaling(kFALSE)) ;

where I have previously included the ratio in the data:

RooFormulaVar ratio_func("ratio","@0/@1",RooArgList(*x,*y));
data->addColumn(ratio_func) ;

But this gives me an Error:

[#0] ERROR:InputArguments -- RooProfileLL::nll_model_data_Profile[ratio]:fillHistogram: WARNING: variable is not an explicit dependent: ratio [#0] ERROR:Plotting -- RooProfileLL::nll_model_data_Profile[ratio]:createPlotProjection: "ratio" is not a dependent and will be ignored.

Is there a way to plot the ratio after the profile likelihood?

Thanks,
Anthony

7

Hi,

I realized only belatedly that there was still a question in your last posting.
The reason that you get the error is that your pdf is (as far as I can tell)
still defined in terms of x and y, even if you have defined a new observable
ratio=x/y in your data.

To be able to make a profile likelihood you need a new pdf that describes
your data in terms of ratio (and either x or y). To do so you can e.g. replace
your observable x in your pdf expression with a function expression ‘ratio*y’,
where ratio is the observable in your dataset.

Wouter

1 Like
8

hi
I need help and I hope that you can help me; I am new with roofit
I am trying to do this



  TString inner  ("(((gamma_CaloHypo_CellID&12288)>>12)== 2)");
  TString middle ("(((gamma_CaloHypo_CellID&12288)>>12)== 1)");
  TString outer  ("(((gamma_CaloHypo_CellID&12288)>>12)== 0)");
  TString conv   ("(gamma_Converted == 1)");
  TString noconv ("(gamma_Converted == 0)");
  TString gfac ("(");
  gfac += "0.976*("+outer+"*"+noconv+")";
  gfac += "+0.994*("+middle+"*"+noconv+")";
  gfac += "+1.001*("+inner+"*"+noconv+")";
  gfac += "+0.959*("+outer+"*"+conv+")";
  gfac += "+0.999*("+middle+"*"+conv+")";
  gfac += "+1.019*("+inner+"*"+conv+")";
  gfac += ")";


  

  RooFormulaVar gFac("gfac1", gfac, RooArgList(gamma_CaloHypo_CellID, gamma_Converted));


  RooDataSet *dataB = new  RooDataSet("dataB","dataB",t,SetOfVariables);


  dataB->addColumn(gFac) ;

    
  dataB->Print("v");

when I do the print what is get is


   39)                      gFac = 0 C  L(-INF - +INF)  ""(0.976*((((gamma_CaloHypo_CellID&12288)>>12)== 0)*(gamma_Converted == 0))+0.994*((((gamma_CaloHypo_CellID&12288)>>12)== 1)*(gamma_Converted == 0))+1.001*((((gamma_CaloHypo_CellID&12288)>>12)== 2)*(gamma_Converted == 0))+0.959*((((gamma_CaloHypo_CellID&12288)>>12)== 0)*(gamma_Converted == 1))+0.999*((((gamma_CaloHypo_CellID&12288)>>12)== 1)*(gamma_Converted == 1))+1.019*((((gamma_CaloHypo_CellID&12288)>>12)== 2)*(gamma_Converted == 1)))""

i.e it adds it as a string

what should I do? what is that I am doing wrong

please help :smiley:

thanks in advance

Cheers