Issue with a simple macro to Fit TH1

[dyba]

Hi guys…
I found this strange behavior while trying to fit a TH1 in a simple macro I wrote. I will try to explain…

I have some TH2Ds. I want to do the ProjectionY on these TH2s, fit the projection with a gauss function and fill a TH1D with the mean obtained from the gauss fits. Finally, I need to fit this TH1D with a pol1 (see image).

In the figure, you can see the results of these operations with the code I wrote. The red triangles are the bin contents (with errors) of the TH1D filled with the mean value from the gauss fit.

As you can see, it basically works, but I don’t know why I have those points at x > 3000 (considering that in my macro there is a check on the entries of the ProjectionY, if they are less than 20 the bin is skipped).

My program should not fill those points and in fact, if I stop the cursor over them, the histogram name is not present.
Then the real problem is when I try to fit the red TH1D with a pol1.
It seems that the fit - range defined between 0 and 2800 - is forced to cross the red triangle at y = 0 and x ~ 2700, and this of course makes the fit wrong.
If I change the fit range excluding the last point, the fit is correct but unfortunately this is not a viable option considering what I need.

I don’t know what are these points at y=0… and I don’t know why the fit seems forced to pass over these points.
I also read that bins with error = 0 should be excluded from a fit, so why do I see this behavior?
How could I fix it?

I hope I was clear in the explanation.
Thank you very much for any help

[Dilicus]

Hi @dyba and welcome to the ROOT forum.
I think filling a TGraph or a TGraphErrors instead of TH1D is better for what are you doing.

In any case a minimal working reproducer of your problem could be useful to understand what is happening.

Best
Stefano

[dyba]

Hi @Dilicus,
thanks for the welcome. I will try with the TGraphErrors. Anyway, this is a small part of my macro with the lines regarding this issue (but the TH2D is in an external root file…)

void do_test(){

   TFile * _file0 = new TFile("root_files/preSel_v02_L7+_CT3_vetoCustom_2018_08-12.root","read");
   TH2D * h_veto;
   TH1D * h_fit_gauss;
   Int_t bin_le,vi,pi,ci;
   TH1D *h_profY = NULL;
   Int_t bin_max;
   TSpectrum s;
   int nfound;
   int nPeaks;
   double *peakPositions;
   double binWidth;
   TString h_name = Form("veto/notHe/veto_vs_LYSO/h_vetoSigHG_P%d-%d_vs_Cry_%d_notHe_LysoHit", 0,0,0);
   h_veto = (TH2D*)_file0->Get(h_name);
   TString h_name3 = Form("h_fit_gauss_%d-%d_vs_Cry_%d", 0, 0, 0);
   h_fit_gauss = new TH1D(h_name3, h_name3,
			  h_veto->GetNbinsX() / 11,
			  0, h_veto->GetXaxis()->GetXmax());
   h_fit_gauss->SetMarkerSize(1.);
   h_fit_gauss->SetMarkerStyle(23);
   h_fit_gauss->SetMarkerColor(2);
   h_fit_gauss->SetLineColor(2); 
   
   Double_t y_max_value_vs_lyso;  
   Double_t x_max_value_vs_lyso;  
   Double_t y_fit_gauss_vs_lyso;  
   Double_t y_fit_gauss_err_vs_lyso; 
   int lyso_fit_max_lateral[9] = {2800, 1700, 2200, 1800, 2900, 2900, 3000, 2700, 2200};
   
   int line = 0;
   for(bin_le=1; bin_le <=  h_veto->GetNbinsX(); bin_le+=11)  {
      h_profY = h_veto->ProjectionY("_temp",bin_le,bin_le+10);
      if(h_profY->GetEntries() < 20)
	 continue;
      bin_max = h_profY->GetMaximumBin();
      nfound = s.Search(h_profY, 2, "", 0.05);
      nPeaks = s.GetNPeaks();
      peakPositions = s.GetPositionX();
      binWidth = h_profY->GetBinWidth(1);
      for (int ip = 0; ip < nPeaks; ip++) 
      {
	 if (peakPositions[ip] > -10 &&
	     peakPositions[ip] <  80   )
	 {
	    h_profY->Fit("gaus","","", peakPositions[ip]-binWidth*6, peakPositions[ip]+binWidth*6);
	    TF1 *f = h_profY->GetFunction("gaus");
	    y_fit_gauss_vs_lyso = f->GetParameter(1);
	    y_fit_gauss_err_vs_lyso = f->GetParError(1);
	    if((peakPositions[ip] - y_fit_gauss_vs_lyso ) < 10){
	       
	       std::cout << "y_fit_gauss_vs_lyso: " << y_fit_gauss_vs_lyso << std::endl;
	       x_max_value_vs_lyso = h_veto->GetXaxis()->GetBinCenter(bin_le+5);
	       h_fit_gauss->Fill(x_max_value_vs_lyso, y_fit_gauss_vs_lyso);
	       h_fit_gauss->SetBinError(bin_le+5, y_fit_gauss_err_vs_lyso);
	    }
	    break;
	 }
      }
   }
   h_fit_gauss->Fit("pol1","", "", 0, lyso_fit_max_lateral[0]);
   TF1 *f1 =  h_fit_gauss->GetFunction("pol1");
   delete h_profY;
   
}

[dyba]

Ok, the TGraphErrors solved the problem. Thanks!

But I don’t understand why the TGraph (in black) and TH1D (in red) have different errors, they are filled with the same content. I just added the lines:

g_fit_gauss->SetPoint(g_fit_gauss->GetN(), x_max_value_vs_lyso, y_fit_gauss_vs_lyso);

g_fit_gauss->SetPointError(g_fit_gauss->GetN()-1, 0, y_fit_gauss_err_vs_lyso);

Still curious to know what is the problem with my first version.

[Dilicus]

It is hard to guess without a file but a I have a couple of suggestions

• I will use a “<” instead of “<=”, I think you do a projection outside the X range of the histogram, and I do not know how ProjectionY works in this case, maybe is the reason for the point with zero error

@dyba

for(bin_le=1; bin_le <= h_veto->GetNbinsX(); bin_le+=11)

• below i will try to use SetBinContent instead of Fill with a weight, just to be consistent. Not 100% sure but maybe this affects how the chi2 is computed.

@dyba

h_fit_gauss->Fill(x_max_value_vs_lyso, y_fit_gauss_vs_lyso);

Best

[dyba]

I think you are right. By replacing the Fill instruction with SetBinContent it seems to work as expected.
Anyway, the TGraphErrors is probably the best option, I don’t know why I didn’t think about it.

Thanks again!