Fitting many histogram and save parameters

diboye · March 13, 2017, 1:52pm

Hi Everyone,
I have many histo that I would like to fit and save their parametters in a txt file and something like that works:

 TFile* f0   = new TFile("nominal_signalDistributions_20170217.root");
    TH1F * chan_4e_15GeV   = (TH1F*)f0->Get("chan_0_hist_150");
    TH1F * chan_4e_20GeV   = (TH1F*)f0->Get("chan_0_hist_200");
    TH1F * chan_4e_25GeV   = (TH1F*)f0->Get("chan_0_hist_250");
    TH1F * chan_4e_30GeV   = (TH1F*)f0->Get("chan_0_hist_300");
    TH1F * chan_4e_35GeV   = (TH1F*)f0->Get("chan_0_hist_350");
    TH1F * chan_4e_40GeV   = (TH1F*)f0->Get("chan_0_hist_400");
    TH1F * chan_4e_45GeV   = (TH1F*)f0->Get("chan_0_hist_450");
    TH1F * chan_4e_50GeV   = (TH1F*)f0->Get("chan_0_hist_500");
    TH1F * chan_4e_55GeV   = (TH1F*)f0->Get("chan_0_hist_550");
    TH1F * chan_4e_60GeV   = (TH1F*)f0->Get("chan_0_hist_600");
    //chan_4e_15GeV->Fit("gaus");
TF1 *fit = new TF1("fit","gaus");
FILE *fp = fopen("filename.dat","w");
chan_4e_15GeV->Fit(fit);
if (fp!=NULL) {
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
       fprintf(fp,"%d %f",i,value);
   }
}

from this I end up to do something that works but not in a proper way since I have to repeat the for loop 10 times which is the number of the histograms.
Could someone show me how to do this in a proper way please.
Cheers

couet · March 13, 2017, 2:08pm

{
   TFile* f0   = new TFile("nominal_signalDistributions_20170217.root");
   TF1 *fit = new TF1("fit","gaus");
   FILE *fp = fopen("filename.dat","w");
   if (fp==NULL) return;
  
   int n = 150;
   TH1F *h;
   for (int i=0; i<10 ; i++) { 
      h = (TH1F*)f0->Get(Form("chan_0_hist_%d",n));
      h->Fit(fit);
      for (int i=0;i<fit->GetNpar();i++) {
         Float_t value = fit->GetParameter(i);
         fprintf(fp,"%d %f",i,value);
      }
      n = n+50;
   }
}

diboye · March 13, 2017, 4:33pm

Hi,
Thanks for replying, your code seems to be what I really want, but I am having some segmentation fault.
Cheers

Axel · March 13, 2017, 8:01pm

Building on top of Olivier’s reply:

// file fitall.C
void fitall() {
   TFile* f0   = new TFile("nominal_signalDistributions_20170217.root");
   if (!f0) {
      std::cerr << "Cannot open file nominal_signalDistributions_20170217.root\n";
      return;
   }

   TF1 *fit = new TF1("fit","gaus");
   FILE *fp = fopen("filename.dat","w");
   if (fp==NULL) return;
  
   for (int n = 150; n <= 600 ; n += 50) { 
      TH1 *h = 0;
      f0->GetObject(Form("chan_0_hist_%d", n), h);
      if (!h) {
         std::cerr << "Cannot get histogram chan_0_hist_" << n << '\n';
         return;
      }
      h->Fit(fit);
      for (int ipar=0; ipar < fit->GetNpar(); ipar++) {
         Float_t value = fit->GetParameter(ipar);
         fprintf(fp,"%d %f", ipar, value);
      }
      fprintf(fp, "\n");
   }
}

That should help you find the bug

Axel.

couet · March 14, 2017, 8:00am

Yes sorry I wrote on the fly without testing … but it was “the idea”…

diboye · March 14, 2017, 12:32pm

Hi Axel,
Thanks for replying and here is what it wrote:

Cannot get histogram chan_0_hist_150
However my code works, it’s here:

void signalFit()
{ TFile* f0   = new TFile("nominal_signalDistributions_20170217.root");
  TFile* f1   = new TFile("FitResult_4e.root","recreate");
    TH1F * chan_4e_15GeV   = (TH1F*)f0->Get("chan_0_hist_150");
    TH1F * chan_4e_20GeV   = (TH1F*)f0->Get("chan_0_hist_200");
    TH1F * chan_4e_25GeV   = (TH1F*)f0->Get("chan_0_hist_250");
    TH1F * chan_4e_30GeV   = (TH1F*)f0->Get("chan_0_hist_300");
    TH1F * chan_4e_35GeV   = (TH1F*)f0->Get("chan_0_hist_350");
    TH1F * chan_4e_40GeV   = (TH1F*)f0->Get("chan_0_hist_400");
    TH1F * chan_4e_45GeV   = (TH1F*)f0->Get("chan_0_hist_450");
    TH1F * chan_4e_50GeV   = (TH1F*)f0->Get("chan_0_hist_500");
    TH1F * chan_4e_55GeV   = (TH1F*)f0->Get("chan_0_hist_550");
    TH1F * chan_4e_60GeV   = (TH1F*)f0->Get("chan_0_hist_600");
    //chan_4e_15GeV->Fit("gaus");
TF1 *fit = new TF1("fit","gaus");
FILE *fp = fopen("Chan_4e.dat","w");
chan_4e_15GeV->Fit(fit);
 chan_4e_15GeV->Write();
if (fp!=NULL) {
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+15,value);//will write for instance 15 1.81947416 14.90458417 0.220540 (15 means 15 GeV)
   }
}
chan_4e_20GeV->Fit(fit);
chan_4e_20GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+20,value);//will write 20 1.87296621 19.88915322 0.291102 (20 means 20 GeV)
   }
   }

chan_4e_25GeV->Fit(fit);
chan_4e_25GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+25,value);
   }
   }
chan_4e_30GeV->Fit(fit);
chan_4e_30GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+30,value);
   }
   }

chan_4e_35GeV->Fit(fit);
chan_4e_35GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+35,value);
   }
   }
chan_4e_40GeV->Fit(fit);
chan_4e_40GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+40,value);
   }
   }

chan_4e_45GeV->Fit(fit);
chan_4e_45GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+45,value);
   }
   }

chan_4e_50GeV->Fit(fit);
chan_4e_50GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+50,value);
   }
   }

chan_4e_55GeV->Fit(fit);
chan_4e_55GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+55,value);
   }
   }
chan_4e_60GeV->Fit(fit);
chan_4e_60GeV->Write();
if (fp!=NULL) {
fprintf(fp,"\n");
   for (int i=0;i<fit->GetNpar();i++) {
      Float_t value = fit->GetParameter(i);
      fprintf(fp,"%d %f",i+60,value);
   }
   }






   fclose(fp);
}

but you can see that’s not an advanced coding.
Cheers

Axel · March 14, 2017, 12:43pm

Hi,

Makes no sense Can you post the ROOT file so we can have a look what’s going on?

Axel.

diboye · March 14, 2017, 12:49pm

Hi,
Here you can find the root file:
/afs/cern.ch/work/d/diboye/public/nominal_signalDistributions_20170217.root
cheers

Axel · March 14, 2017, 1:00pm

Thanks. Your hists are TH1D not TH1F! I have updated my code above to reflect that - as you don’t care about the actual precision in the code we can just use a pointer to the base class: TH1.

(And claiming it’s a TH1F even though it’s a TH1D can actually give wrong results, so it’s good we found out…)

Cheers, Axel.

couet · March 14, 2017, 1:00pm

when I run this macro with the root file you just sent,

{
   TFile* f0   = new TFile("nominal_signalDistributions_20170217.root");
   TF1 *fit = new TF1("fit","gaus");
   FILE *fp = fopen("filename.dat","w");
   if (fp==NULL) return;

   int n = 150;
   TH1F *h;
   for (int i=0; i<10 ; i++) {
      h = (TH1F*)f0->Get(Form("chan_0_hist_%d",n));
      h->Fit(fit);
      for (int i=0;i<fit->GetNpar();i++) {
         Float_t value = fit->GetParameter(i);
         fprintf(fp,"%d %f\n",i,value);
      }
      n = n+50;
   }
}

… I get the following file:

0 0.584926
1 14.700195
2 0.426157
0 0.636955
1 19.615734
2 0.538808
0 0.604746
1 24.563480
2 0.669224
0 0.639619
1 29.402952
2 0.807544
0 0.661607
1 34.252602
2 0.926711
0 0.791310
1 39.253941
2 1.030971
0 0.789471
1 44.173283
2 1.180475
0 0.787339
1 49.006859
2 1.363823
0 0.865696
1 53.964760
2 1.504936
0 1.000655
1 58.901871
2 1.415896

diboye · March 14, 2017, 1:34pm

Hi Guys,
Thanks a lot your codes looks really impressing comparing to mine.
With Axel’s code I got the same the result however I have not yet end up the reproduce Couet’s result we could try to understand why for someone else who will have a such problem; in fact I got this:
Error: Symbol lt is not defined in current scope fittingRtk.C:46:
Error: Binary operator oprand missing fittingRtk.C:46:
I am puting the entire code here:

#include "TROOT.h"    
#include "TObject.h" 
#include "TFile.h"
#include "TTree.h"
#include "TH2.h"
#include "TMath.h"
#include "TGraph.h"
#include "TColor.h"
#include "TCanvas.h"
#include "TLegend.h"
#include <TStyle.h>
#include <iostream>
#include "TLatex.h"
void fitall() {
 TFile* f0   = new TFile("nominal_signalDistributions_20170217.root");
   TF1 *fit = new TF1("fit","gaus");
   FILE *fp = fopen("filename.dat","w");
   if (fp==NULL) return;
  
   int n = 150;
   TH1F *h;
   for (int i=0; i&lt;10 ; i++) { 
      h = (TH1F*)f0-&gt;Get(Form("chan_0_hist_%d",n));
      h-&gt;Fit(fit);
      for (int i=0;i&lt;fit-&gt;GetNpar();i++) {
         Float_t value = fit-&gt;GetParameter(i);
         fprintf(fp,"%d %f \n",i,value);
      }
      n = n+50;
   }

}

cheers

couet · March 14, 2017, 1:44pm

Let me attached the .C file. The message you get shows that you did a wrong copy/paste. The symbols > and < have been replaced … look at your code it is full of &lt and &gt …

diboye.C (469 Bytes)

system · March 28, 2017, 1:44pm

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