Two types of an user-defined fit function

Sisisan · September 11, 2017, 3:29pm

Hi All,

Could you help me, please, to comprehend what the differences are the following formulas for fiting ?

First is here :

Double_t fitfunc(Double_t *x, Double_t *par) 
{
Double_t f;
f=par[0]*TMath::Power(x[0],(par[1]*10));
return f;
}

TF1 *fit = new TF1("fit",fitfunc,550,800,2); // here, access parameter number 
fit->SetParameters(0.000001,0.75);
fit->SetParNames("Constant", "k value");

Second is here :


TF1 *fit = new TF1("fit","[0]*TMath::Power(x,([1]*10))",550,800); // here, don't access parameter number
fit->SetParameters(0.000001,0.75);
fit->SetParNames("Constant", "k value");

Cheers,
Ersel

Wile_E_Coyote · September 11, 2017, 3:31pm

Double_t fitfunc(Double_t *x, Double_t *par)

f = par[0] * TMath::Power(x[0], (par[1] * 10.));

fit->SetParNames(“Constant”, “k value”);

Sisisan · September 11, 2017, 3:38pm

Ok, thank you. I wrote those accidentally. Now I fixed them. Anything else ?

Cheers,
Ersel

Wile_E_Coyote · September 11, 2017, 3:39pm

You still haven’t fixed all bugs that I pointed out.
Compare my lines with yours, character by character.

Sisisan · September 11, 2017, 3:43pm

OK ,done, i think

Wile_E_Coyote · September 11, 2017, 3:44pm

So, now, after you fixed it … what’s your problem?

Sisisan · September 11, 2017, 4:08pm

I have many problems for fitting with this (above) formula ,which is a theoretical formula for fitting to my data and is correct i’m sure %100

First problem is here :

What does “Invalid FitResult” mean ? Does it come out due to Chi2/Ndf >> 1 ?

p.s. : while some fit results as in the following have “invalid result”, some of them have not !

****************************************
         Invalid FitResult  (status = 4 )
****************************************
Minimizer is Minuit / Migrad
Chi2                      =      635.553
NDf                       =            2
Edm                       =       38.969
NCalls                    =         1224
Constant                  =  2.58789e-15   +/-   1.60865e-15
k value                   =      0.71379   +/-   0.00880427
0.71379

Second is here : Why does “k value” in the formula ,which is slope, change always as the starting points change ?

For instance,

SetParameters(0.000000000000001,0.8); k value of fit results is 0.847833 as in the following results

**************************************
Minimizer is Minuit / Migrad
Chi2                      =      14.0821
NDf                       =            2
Edm                       =  2.74101e-09
NCalls                    =          863
Constant                  =  3.81246e-19   +/-   8.2229e-19
k value                   =     0.847833   +/-   0.0327148

SetParameters(0.000000000000001,0.7); k value of fit results is 0.846402 as in the following results

****************************************
         Invalid FitResult  (status = 4 )
****************************************
Minimizer is Minuit / Migrad
Chi2                      =      14.2095
NDf                       =            2
Edm                       =      3.07366
NCalls                    =         1223
Constant                  =  4.18651e-19   +/-   1.24307e-19
k value                   =     0.846402   +/-   0.00427859

Also, i have some problems for this issue. I’ll post them here as soon as i meet them

Cheers,
Ersel

Wile_E_Coyote · September 11, 2017, 4:32pm

I think you need to attach your histogram (or graph) which your try to fit.

Sisisan · September 11, 2017, 5:01pm

Fitted histogram and result is here

Minimizer is Minuit / Migrad
Chi2 = 14.0821
NDf = 2
Edm = 2.98678e-07
NCalls = 765
Constant = 3.81232e-19 +/- 7.98971e-19
k value = 0.847833 +/- 0.0317883

Uploading…

Wile_E_Coyote · September 11, 2017, 5:14pm

I think the ROOT file with the original histogram is needed (plus a small macro which fits it).

Sisisan · September 11, 2017, 5:31pm

Macro is here , which calculates GAIN of a PMT and draws GAIN vs HV, and then fits GAIN vs HV

to run : root -l ‘BURNINFitforOnePMT.C(2949)’

using namespace std;
char rootfile[200];

Double_t func(Double_t *x, Double_t *par) {

  Int_t n = 10;
  Double_t f;
  f = par[0]*TMath::Power(x[0],(par[1]*n));
  return f;

}

void BURNINFitforOnePMT (Int_t runno) {

  Double_t Constant;
  Double_t k;
  Double_t k_mean;

  TH1D *kdist = new TH1D("kdist","k distribution",100,0.5,0.9);
  TCanvas *c2 = new TCanvas("c2","k distribution",600,800);
  TCanvas *c1 = new TCanvas("c1","Gain vs HV",600,800);
  TH2D *GAINHV  = new TH2D("h1","GAINvsHV",1000,550,800,1000,0,3000000);

  TF1 *fit = new TF1("fit1",func,550,800,2);

  sprintf(rootfile,"/Users/core/Desktop/!!!THESIS!!!/PMT ANALYSIS /BURNIN Thesis/BURNINat4HVs/LED_4HVs_ROOTFILES/Run_%d.root",runno);
  TFile *f = new TFile(rootfile);

 TFile *f1 = new TFile("GAINHVFIT.root","RECREATE");

  ifstream in;
  in.open("/Users/core/Desktop/!!!THESIS!!!/PMT ANALYSIS /BURNIN Thesis/txt/B904/map.txt",ios::in);

  Int_t HV[] = {600,650,700,750};
  Int_t qie,ieta,depth;

  while(in>>qie>>ieta>>depth) {

     if(qie==0001 && ieta==33 && depth==1) {

      for(Int_t db=0; db<1; db++) {

        for(Int_t hv=0; hv<4; hv++) {

        stringstream pedname;
        stringstream ledname;
        TString HFLEDname;
        TString HFPEDname;
        HFLEDname.Form("/Analyzer/QIEsSumLED%d_QiECh_%d_DBOX_%d_eta_%d_D_%d",HV[hv],qie,db,ieta,depth);
        HFPEDname.Form("/Analyzer/QIEsSumPED%d_QiECh_%d_DBOX_%d_eta_%d_D_%d",HV[hv],qie,db,ieta,depth);
        pedname<<HFPEDname;
        ledname<<HFLEDname;

        TH1F *histled = (TH1F*)f->Get(ledname.str().c_str());
        TH1F *histped = (TH1F*)f->Get(pedname.str().c_str());

        ledname.str("");
        pedname.str("");

        if(histled && histped)  {

          Double_t LEDMEAN;
          Double_t PEDMEAN;
          Double_t LEDRMS;
          Double_t PEDRMS;

          LEDMEAN = histled->GetMean();
          PEDMEAN = histped->GetMean();
          LEDRMS = histled->GetRMS();
          PEDRMS = histped->GetRMS();

          Double_t TRUEMEAN = LEDMEAN - PEDMEAN;
          Double_t TRUERMS = sqrt(pow(LEDRMS,2)-pow(PEDRMS,2));
          Double_t NPE = 1.15*(pow(TRUEMEAN,2)/pow(TRUERMS,2));// 1.15+-0.05 = ENF(Excess Noise Factor)
 //Double_t NPE = 1.15*(pow(LEDMEAN,2)/pow(LEDRMS,2));
          Double_t Charge = NPE*1.6*pow(10,(-19));
          Double_t GAIN = (TRUEMEAN*2.6)*pow(10,(-15))/Charge;
          //Double_t GAIN = LEDMEAN/NPE;
          //Double_t GAIN = LEDMEAN;


            if(GAIN>0) {


              GAINHV->Fill(HV[hv],GAIN);
              GAINHV->SetMarkerColor(2);
              GAINHV->SetMarkerStyle(8);

              fit->SetParameters(0.000000000000001,0.9);
              fit->SetParNames("Constant","k value");

            }
        }
        }
        }
 c1->cd();
      GAINHV->Fit("fit1");
      //      GAINHV->Draw();


      TF1 *fitResults = GAINHV->GetFunction("fit1");
      Constant = fitResults->GetParameter(0);
      k = fitResults->GetParameter(1);

      kdist->Fill(k);

      }
  }

  c2->cd();
  kdist->Draw();

  k_mean = kdist->GetMean();
  cout<<k_mean<<endl;

  c2->Write();
  c1->Write();

}

And rootfile is in the attachment.

And I also notice that the fit result changes as ‘histogram bin number’ of TH2D changes ,which is 1000 in the script !?

Run_2949.root (774.1 KB)

Cheers,
Ersel

Wile_E_Coyote · September 11, 2017, 5:42pm

map.txt file is missing.

Sisisan · September 11, 2017, 5:46pm

oops ok,sorry, it’s also here ;
map.txt (241 Bytes)

Wile_E_Coyote · September 11, 2017, 7:23pm

Try with:
f = TMath::Power((par[0] * x[0]), (par[1] * n));
and:
fit->SetParameters(7.e-3, 0.8);
or with:
f = TMath::Power((x[0] / par[0]), (par[1] * n));
and:
fit->SetParameters(150., 0.8);

Sisisan · September 12, 2017, 6:37pm

It works very well, thank you very much

system · September 26, 2017, 6:37pm

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