Fit NKG distribution

Hi, I wrote the macro in the attachment, it must draw a plot and fit the histograms using the NKG distribution

where r0=100 is the moliere radius; it compiles, but the fits aren’t good, because I’ve

but I must have something like this

What is the problem?
distrlat.c (5.39 KB)


   TF1 *fitdistrlat = new TF1("fitdistrlat", "[0]/10000.*TMath::Power((x/100.), ([1]-2.))*TMath::Power((1.+x/100.), ([1]-4.5))", 1, 1e4);
   fitdistrlat->SetParameters(1e-4, -8e-1);

Hi pepe thanks, ufortunatly this is the result:

Befeore looking your reply I changed the code as in the attachment, but this was the result:

maybe it can be a problem setting the parameters (I don’t know their values…)
distrlat.c (6.14 KB)


   TF1 *fitdistrlat = new TF1("fitdistrlat", "[0]/[2]/[2]*TMath::Power((x/[2]), ([1]-2.))*TMath::Power((1.+x/[2]), ([1]-4.5))", 1, 1e4);
   fitdistrlat->SetParNames("cN", "s", "r0");
   fitdistrlat->SetParameters(1., 1., 100.);
   // fitdistrlat->FixParameter(2, 100.); // fix "r0"

Hi pepe…no… unfortunatly it doesn’t draw the fit lines like in this plot

PS. This is the .root file … .root?dl=0

You can try the attached macro, but it seems to me that your data do not really like your function.
distrlat.cxx (5.74 KB)

Hi pepe…I know that usually the lateral distribution of an EAS is described by the NKG function, so I tried to fit the data using this function (but the supervisor is on holiday, so I didn’t ask him if it is correct)…
I runned your macro and this is the result:

it looks like better than other macros…
what do you think about?

I will try the macro for other set of data too (I’ve 10 CORSIKA file, 5 files are EAS induced by photons and 5 files by protons)

I tried the code to others data file…but for some files the results are worse…maybe as you wrote it’s a problem of data…

Maybe your “analytical function” is not really describing your “experimental data” well (you would need to ask corsikans what function they are actually using).

Try also:

  TFile *f = TFile::Open("cors_plot-gamma1.root");
  TH1F *hgamlat = (TH1F *)f->Get("hgamlat");
  TH1F *helelat = (TH1F *)f->Get("helelat");
  TH1F *hmuolat = (TH1F *)f->Get("hmuolat");
  TH1F *hhadlat = (TH1F *)f->Get("hhadlat");
  TCanvas *c = new TCanvas("c", "c");
  c->Divide(2, 2);
  c->cd(1); gPad->SetLogx(); hgamlat->Draw("E1");
  c->cd(2); gPad->SetLogx(); helelat->Draw("E1");
  c->cd(3); hmuolat->Draw("E1");
  c->cd(4); hhadlat->Draw("E1");

Note that, for all your histograms, “bin_error” is something like 100 to 100000 times bigger than “bin_content” (I would expect the ratio “bin_error / bin_content < 1”).
Is it possible that these histograms were scaled, but Sumw2 was not called before making this operation?

Attached is a “cleaned” version of your macro (fits the “NKG” function).
It should be easier to modify / (re)use it.
A variation of this macro, which uses the “Linsley” function instead of “NKG”, can be found here.
distrlat.cxx (5.17 KB)

Hi pepe, i tried the macros, but with same results…maybe you are right that the data are not described by that function

Actually, I started to think that the main problem may be related to “bin errors” of all your histograms (which I pointed out in my previous post). Then the chi^2 is ridiculously small, which may fool the minimizer.
BTW. The current macro uses the same (fixed) “r_mol” parameter for all types of particles and I guess it’s not really completely right (but it depends on what was used by the one who produced this simulation). One could also try to fit “r_mol” (one simply needs to comment out the “FixParameter” line in the macro).

Hi pepe, I wrote to the author of simulations to know the function used for the simulations, but he didn’t reply…