Fitting histograms with Roofit inside "for loop"CODE CRASH!

Roofit and RooStats
1

Hi,
I am a new programmer and trying to use RooFit for the first time to fit a landau curve on a histogram. I have used MakeClass to create a .h and a .C file from my original root file. Hence I am using the macro in the .C file for filling my histogram. Initially I chose an individual histogram and imported the data of the histogram using the code:

    //importing the histogram data
     RooRealVar x("x","x",-0.5,100.5);
TH1F *hh=(TH1F*)gDirectory->Get("set2"); //set2 is the histogram name
RooDataHist *dh=new RooDataHist("dh","data histogram",x,hh);

    //creating a landau pdf
    RooRealVar ml("mean","mean of landau",x1);
    RooRealVar sl("sigma","sigma of landau",0.5,0.445,0.523);
    RooLandau landau("landau","landau pdf",x,ml,sl);

   RooPlot *frame1=x.frame(Title("adc1_vs_count_set2_1757"),Bins(101));
   dh->plotOn(frame1);
   landau.fitTo(*dh,Range(40.5,70.5));
   landau.plotOn(frame1);

TCanvas *c_set2=new TCanvas(“c_set2”,“adc1_vs_count”,1000,1000);
c_set2->SetGrid();
frame1->SetXTitle(“count”);
frame1->SetYTitle(“adc1_raw”);
frame1->Draw();

This code worked nicely and the histogram was properly imported and the landau pdf was fitted nicely.

But now I want to fit 20 histograms with the landau pdf and would like to run a for loop for this.
I tried using the code as follows:

ostringstream  var_set2,varml_set2,varsl_set2,varlan_set2 varframe_set2;
    RooDataHist *dh[40]
{
TCanvas *c_set2=new TCanvas("c_set2","adc_count",1000,1000);
c_set2->Divide(5,4);
    
    for(int ievt=1;ievt<21;ievt++)

{

  RooRealVar x("x","x",-0.5,100.5);
   
		
 vardh_set2<<"dh"<<ievt;
     dh[ievt-1]=new RooDataHist(vardh_set2.str().c_str(),"data hist",x,h2b_set2[ievt-1]);//!!!h2b_set2[40] are the array of histograms declared before
      vardh_set2.str("");
	
	
	RooRealVar ml("ml","landau mean",53.0);
	RooRealVar sl("sl","landau sigma",1.0);
	RooLandau landau("landau","landau pdf",x,ml,sl);
	
	
	
       varframe_set2<<"xframes"<<ievt;
	RooPlot *xframe=x.frame(Name(varframe_set2.str().c_str()),Title("framework"),Bins(101));
	varframe_set2.str("");
	
	
	dh[ievt-1]->plotOn(xframe);//!!!!!!!!!!THIS IS WHERE MY CODE CRASHES DURING THE FIRST ENTRY!
	landau.fitTo(*dh[ievt-1],Range(40.5,70.5));
           landau.plotOn(frame1);
	
          c_set2->cd(ievt);
	
	
	xframe->Draw();
	c_set2->Update();
}

When I tried running this the code crashed. I might have made serious mistakes in the code as I am very naive with Root. Basically I wanted to know how to import 20 histograms and create plot frames for each individual histograms in a “for loop” and fit them with the same landau pdf and draw them on different pads on a canvas using RooFit. Eventually I have to fit 90000 adc histograms. So I just tried to run a for loop for 20 histograms first. Please help.
timewalk2bset2.C (7.86 KB)

2

Hi,

I think you need to have the variable “x” associated too the data set declared outside the loop. Otherwise it will be deleted after the first iteration and the first data set will have lost the link to the corresponding variable.
Alternatively you can create 20 different variables.

Please post next time also your data file, so one can try to run your code

Best Regards

Lorenzo

3

Thank you for your response. The data file is 1.3 GB so I was unable to attach the file with my previous post. I will try the solution that you suggested and try to run the code.

Regards,
Christine

4

Hi,
I am still not able to create separate RooRealVar variables for individual histograms in a for loop as was suggested in the last post. And also how do I create separate plot frames for each histogram inside the loop. I have maked in the code where I have made changes. I still could not attach the data file since its too huge. Please help me with the code.
timewalk2bset2.h (5.72 KB)
timewalk2bset2.C (6.6 KB)

5

Hi,
Cannot you not post the file somewhere, like /afs or dropbox or whatever that can be downloaded ?
Or eventually you could reduce it to have a smaller number of events, or divide the macro in two: first make the histograms from the tree, save it in a file and later you have a macro doing the fitting. Then you send me the file with the histograms and the macro doing the fitting.
Without running the code it is difficult to find the problem.

Best Regards

Lorenzo

6

Hi Lorenzo,

I can upload the root data file as well as the other relevant files in dropbox or google drive. Please let know how can I share the files with you.

Regards,
Christine

7

Hi Christine,
Please update to dropbox and you can send me by private message (at my cern e-mail address) the links to the file

Regards

Lorenzo

8

Hi,

Thank you for the ROOT file. How do you run exactly the code and do the crash ?
It seems to me that t…::Loop() does not loop on all the entries in the tree

Lorenzo

9

Hi Lorenzo,

I was planning to fit the first 20 histograms with the Landau PDF using RooFit. I did that to cross check if I can run RooFit in a for loop. Whne I pass 21 as the arguement to the loop function and run the code, I get a segmentation fault. Not all events has a signal. In the first 20 events only event number 18 has a signal. So if I do the plot for only the event number 18, I can see the landau pdf superimposed on the histogram, but then I think the fit has normalization issues. I am not sure how do I normalize the histogram to get a correct fit.

Thanks,
Christine

10

Hi,

Even, t.Loop(1) crashes for me. Can you please tell me how do you make a plot for event 18 ?
Is the macro you have posted correct ? The code looks a bit messy, lots of things are commented…
What I get is :

I am executed 1
I have been executed
I am near canvas
i am at roodatahist1
[#1] INFO:DataHandling -- RooDataHist::adjustBinning(dh1): fit range of variable x_set21 expanded to nearest bin boundaries: [-0.5,100.5] --> [-0.5,100.5]
i am at rooplot1
I am at dh plotting1
[#0] WARNING:Plotting -- RooHist::addBin(dh1_plot__x_set21) WARNING: negative entry set to zero when Poisson error bars are requested
[#0] WARNING:Plotting -- RooHist::addBin(dh1_plot__x_set21) WARNING: negative entry set to zero when Poisson error bars are requested

 *** Break *** segmentation violation

It seems to me there is a problem in the definition of your bins…

Lorenzo

11

Hi,
There are lot of diffrerent things that I wanted to try in the code so this is the reason why severel things has been commented out. There are 3 for loops in my code. The first for loop declares the histograms and it also declares polynomial zero fits for all the histogram. This polynomial zero fit has been used for offset subtraction since originally all the histograms were starting at 200th mark on the y-axis so this offset subtraction just brings down the histogram to 0 mark. The second for loop is the event loop and here is where I am filling all the histograms. The third for loop is used for fitting the histograms. In the third for loop I have tried to fit the histograms using roofit.

When I try to plot a histogram for a particular event like in this case event number 18.I would make modifications to the code. I would just declare a single histogram and then Fill the histogram inside the event loop with an if statement like
if(eventn_iev==18)
histogram-Fill()
and then fit the histogram using roofit.

I know the code looks a little messy, but since I am very new to programming so I tend to make mistakes a lot.

12

Hi,

Thanks for the explanations. I understand now the code, but not what you really want to do.
The histograms looks to be some time series of raw ADC counts, which can be in some case negatives.
They are not a frequency distribution, thus an empirical estimate of some distribution.
RooFit is used to fit a distribution to a representing data set, you cannot have negative counts and also the pdf cannot be negative !
The crash is probably caused by plotting a data set with negative contents

Lorenzo

13

Hi,

yes the histograms are flash ADC data, which is basically voltage versus time. I eventually have to calculate the rise time of the these adc signals for further analysis. When I am fitting the histogram using roofit for a particular event that has a signal then the code works, I mean when I do it for individual histograms. But as I said that not all events has a signal. Infact the very first histogram for event=1 does have some negative values in the plot. So you are saying that because of the negative contents in the very first event the code crashes? So if I try the roofit for only those events that has a proper signal and fit those in a for loop, then my code would probably work?

14

Yes, I would do some clean up and select first the interesting events before fitting and plotting them.

Lorenzo

15

Hi,
while using roofit, how do I get the maximum y and x values from the fit. Or x values corresponding to a given y value. I mean is there something similar in roofit like TF1::GetX() and TF1::GetMaximum() where we can get the x and the y value from the fit.

Thanks,
Christine

16

What do you mean exactly, how to get the x value corresponding to the maximum of a pdf function obtained from a fit ?

You can always transform a RooFit pdf (a RooAbsPdf) to a TF1 object and the call TF1::GetMaximumX() using the RooAbsPdf::asTF method. Here is a simple example

RooWorkspace w; 
RooAbsPdf * g = w.factory("Gaussian:g(x[-3,3],m[0],s[1])");
TF1 * f1 = g->asTF(*w.var("x"));
cout << f1->GetMaximumX() << endl;

Best Regards

Lorenzo