[Solved] Fit TH1D with Poisson Function

Hi,

I’m trying to fit histogram with Poisson function without any success, where I wrong?

int tempo(){
	
	ifstream file1("muon_Cx.txt");
	
	string nome,stringa,buttare;
	
	double dummy;
 	int conteggio=0;
	
	vector<double> lifetime;
	vector<string> nomi;

	while(getline(file1,stringa)){
	
		nomi.push_back(stringa);
	
		conteggio++;
	
	}
	
	int bin=sqrt(conteggio);

	TH1* h3 = new TH1D("Note", "Distrubuzione dei muoni da raggi cosmici", bin, 0, 10e-6);
	
	for(int i=0;i<conteggio;i++){
		
		vector<double> x,y;
	
		ifstream file2(nomi[i]);
		cout<<nomi[i]<<endl;
		if(file2){cout<<"bene"<<endl;

	}

		string nome=nomi[i];
			
	    for(Int_t i=0;i<6;i++){

	    		getline(file2,buttare);
				
	    	}
	
	   	while(file2>>dummy){
		
	   			x.push_back(dummy);

	   			file2>>dummy;
	   			dummy=-dummy;
	   	
	   			y.push_back(dummy);
	 
	   		}
			
			
   		 TH1* h2 = new TH1D("h2", "h_2", 126, x[0], x[x.size()-1]);
	
   	     for(Int_t i=0;i<x.size();i++){

   		 	h2->Fill(x[i],y[i]);

   		}
		
		
   		TSpectrum *s = new TSpectrum();
		
   		Int_t nfound = s->Search(h2,1);
   		Double_t *xpeaks = s->GetPositionX();
   				
		cout<<abs(xpeaks[1]-xpeaks[0])<<endl;

		h3->Fill(abs(xpeaks[1]-xpeaks[0]));
				
		}
	
		TCanvas *c1=new TCanvas();
	
		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Rate");		
		h3->Draw();

		c1->SaveAs("tempo.pdf");

		TCanvas *c3=new TCanvas();

		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Rate");
        TF1* fit = new TF1("fit","TMath::Poisson(x,[0])",0,10);
        fit->SetParameters(h3->GetMaximum(),h3->GetMean());		
		h3->Fit("fit","R");
		h3->Draw();
		c3->SaveAs("tempo_fit.pdf");		
	
		TCanvas *c2=new TCanvas();

		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Rate");		
		h3->Sumw2();
		h3->Draw();

		c2->SaveAs("tempo_err.pdf");		
	
	return 0;
}

_ROOT Version: 6.13/03
_Platform: Xubuntu 17.10
_Compiler: gcc 7.2.0


You maybe want: "[0] * TMath::Poisson(x, [1])"

BTW. When you post “source code” or “output” here, do remember to enclose them into two lines which contain just three characters ``` (see how your post has been edited above).

c1_n3.pdf (15.2 KB)

Doesn’t fit… i dont understand why…

Try to fit without “R” (maybe the range of your function is not right for your histogram).

Same result… , nothing change.

I think the problem is in the “x range” of your histogram. From your picture, I can see that “x_max = 10.e-6” and the Poisson distribution expects “x_max >= 1” (so, in principle you should also have “mean >= 0.5”).
You could try to “scale” your “x range”:

{
  double xmin = 0.;
  double xmax = 10.e-6;
  TF1 *f = new TF1("f", "[0] * TMath::Poisson((x/[2]), ([1]/[2]))", xmin, xmax);
  double constant = 160.;
  double mean = 1.5e-6;
  double xscaling = 2.e-7;
  f->SetParameters(constant, mean, xscaling);
  f->SetParNames("Constant", "Mean", "XScaling");
  f->SetNpx(1000);
  f->Draw();
  double integral = f->Integral(xmin, xmax);
  std::cout << "integral = " << integral << std::endl;
  std::cout << "scaled integral (constant) = " << integral / xscaling << std::endl;
}

Your code work… but i have some doubts

  1. The ouput is:
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  area         5.17908e+02   8.10529e+01   1.42627e-01  -1.04046e-06
   2  mean         2.18214e-06   1.06192e-07   2.97350e-10   8.55159e+02
   3  scaling      2.41319e+06   3.72806e+05   3.73618e+01   1.84868e-09

How we can see mean life of muon : 2.18214e-06 +- 1.06192e-07 , ok it is into the litterature range, but I don’t understant with which criteria i fix the 3 parameters.

  1. Maybe the fit is wrong, why the maximum of fit is not 160 ?

I re-post the edit code:

#include<iostream>
#include<fstream>
#include<vector>
#include<string>
#include<vector>
#include<math.h>

using namespace std;

int tempo(){
	
	ifstream file1("muon_Cx.txt");
	
	string nome,stringa,buttare;
	
	double dummy;
 	int conteggio=0;
	
	vector<double> lifetime;
	vector<string> nomi;

	while(getline(file1,stringa)){
	
		nomi.push_back(stringa);
	
		conteggio++;
	
	}
	
	int bin=sqrt(conteggio);

	TH1* h3 = new TH1D("Note", "Distrubuzione dei muoni da raggi cosmici", bin, 0, 10e-6);
	
	for(int i=0;i<conteggio;i++){
		
		vector<double> x,y;
	
		ifstream file2(nomi[i]);
		cout<<nomi[i]<<endl;
		if(file2){cout<<"bene"<<endl;

	}

		string nome=nomi[i];
			
	    for(Int_t i=0;i<6;i++){

	    		getline(file2,buttare);
				
	    	}
	
	   	while(file2>>dummy){
		
	   			x.push_back(dummy);

	   			file2>>dummy;
	   			dummy=-dummy;
	   	
	   			y.push_back(dummy);
	 
	   		}
			
			
   		 TH1* h2 = new TH1D("h2", "h_2", 126, x[0], x[x.size()-1]);
	
   	     for(Int_t i=0;i<x.size();i++){

   		 	h2->Fill(x[i],y[i]);

   		}
		
		
   		TSpectrum *s = new TSpectrum();
		
   		Int_t nfound = s->Search(h2,1);
   		Double_t *xpeaks = s->GetPositionX();
   				
		cout<<abs(xpeaks[1]-xpeaks[0])<<endl;

		h3->Fill(abs(xpeaks[1]-xpeaks[0]));
				
		}
	
		TCanvas *c1=new TCanvas();
	
		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Rate");		
		h3->Draw();

		c1->SaveAs("tempo.pdf");

		TCanvas *c3=new TCanvas();

		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Rate");
        TF1 *f = new TF1("f", "[0] * TMath::Poisson([2] * x, [2] * [1])", 0,10e-6);
        double area = 160.;
  		double mean = 1.5e-6;
  		double scaling = 5.e6;
  		f->SetParameters(area, mean, scaling);
  		f->SetParNames("area", "mean", "scaling");
		h3->Fit("f");
		h3->Draw();
		c3->SaveAs("tempo_fit.pdf");		
	
		TCanvas *c2=new TCanvas();

		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Rate");		
		h3->Sumw2();
		h3->Draw();

		c2->SaveAs("tempo_err.pdf");		
	
	return 0;
}

the output is : c1_n25.pdf (15.4 KB)

You can try to play with fit options, e.g.: h3->Fit(f, "L"); // "L" or "WL"

Nothing change. For me the problem is in the parameters…

        double area = 160.;
 		double mean = 1.5e-6;
 		double scaling = 5.e6;
 		f->SetParameters(area, mean, scaling);
 		f->SetParNames("area", "mean", "scaling");

Why do you choose this specific parameter?

I tried to play changing the parameters but the result was the same.

These were just “random” values chosen so that the function’s plot looked more or less similar to your histogram.

Do NOT expect that the “area” parameter is the “height at maximum”, as the Poisson distribution is NOT 1.0 at its “maximum” (and you multiply them).

Sorry but i don’t undertand… I can’t put random numbers in order to obtain a good fit.
For example, can I use member function GetMean() to obtain value of mean? If this is correct which are the other member functions that i need to use?

Then I don’t understand why i need to use 3 parameters and not one, poisson distribution depend by one parameter, the variance.

In ROOT, before you try to fit your graph or histogram, you MUST set “reasonable” initial values for ALL parameters of your function (except for some “built-in” formulas, like “gaus”, for which the standard fit procedure can automatically “guess” them), otherwise the fitting procedure may easily misbehave.

Looking at the result of your fit (in one of your previous posts), I guess this could be fine:

TF1 *f = new TF1("f", "[0] * TMath::Poisson((x/[2]), ([1]/[2]))",
                 h3->GetXaxis()->GetXmin(), h3->GetXaxis()->GetXmax());
// maximum( Poisson( expected value = variance = 5. ) ) = 0.18
f->SetParameters(h3->GetBinContent(h3->GetMaximumBin()) / 0.18, // "Constant"
                 h3->GetMean(), // "Mean"
                 h3->GetMean() / 5.); // "XScaling"
f->SetParNames("Constant", "Mean", "XScaling");
f->SetNpx(1000);
h3->Fit(f, ""); // "" or "L" or "WL" or ...

Now i understand the criteria, but whatever parameters I insert the result does not change, it’s so strange.

Nothing is change, same plot same output…

I understand the error. I make a mistake,is not correct fit my histogramm with Poisson function, because the physics under decay muon can be fitted with an exponential, so i edit my code in this way

#include<iostream>
#include<fstream>
#include<vector>
#include<string>
#include<vector>
#include<math.h>

using namespace std;

int fit_muon(){
	
	ifstream file1("muon_Cx.txt");
	
	string nome,stringa,buttare;
	
	double dummy;
 	int conteggio=0;
	
	vector<double> lifetime;
	vector<string> nomi;

	while(getline(file1,stringa)){
	
		nomi.push_back(stringa);
	
		conteggio++;
	
	}
	
	int bin=sqrt(conteggio);

	TH1* h3 = new TH1D("Note","Distribuzione decadimento muoni cosmici", bin, 0,10e-6);
	
	for(int i=0;i<conteggio;i++){
		
		vector<double> x,y;
	
		ifstream file2(nomi[i]);
		cout<<nomi[i]<<endl;
		if(file2){cout<<"bene"<<endl;

	}

		string nome=nomi[i];
			
	    for(Int_t i=0;i<6;i++){

	    		getline(file2,buttare);
				
	    	}
	
	   	while(file2>>dummy){
		
	   			x.push_back(dummy);

	   			file2>>dummy;
	   			dummy=-dummy;
	   	
	   			y.push_back(dummy);
	 
	   		}
			
				
   		TH1* h2 = new TH1D("h2", "h_2", 126, x[0], x[x.size()-1]);
	
   	    for(Int_t i=0;i<x.size();i++){

   		 	h2->Fill(x[i],y[i]);

   		}
		
		
   		TSpectrum *s = new TSpectrum();
		
   		Int_t nfound = s->Search(h2,1);
   		Double_t *xpeaks = s->GetPositionX();
   				
		cout<<abs(xpeaks[1]-xpeaks[0])<<endl;

		h3->Fill(abs(xpeaks[1]-xpeaks[0]));
				
		}
	
		TCanvas *c3=new TCanvas();

		h3->GetXaxis()->SetTitle("Tempo [ns]");
		h3->GetYaxis()->SetTitle("Conteggio");
		h3->Sumw2();


        TF1 *f = new TF1("f","[0]*exp(-x/[1])",0,10e-6);

		f->SetParameters(160,h3->GetMean());
		f->SetParNames("costant", "meanlife");
		
		h3->Fit(f);
		
		c3->SaveAs("tempo_fit.pdf");		
		
	return 0;
}

It’s works.

Note for future reader : if you want to fit an histogram with Poisson function, this see code is correct.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.