Fit NKG distribution

ROOT
1

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?
Thanks
distrlat.c (5.39 KB)

2

Try:

   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);
   fitdistrlat->SetParName(0,"cN");
   fitdistrlat->SetParName(1,"s");
3

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)

4

Try:

   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"
5

Hi pepe…no… unfortunatly it doesn’t draw the fit lines like in this plot i66.tinypic.com/2wm2nhe.png …

PS. This is the .root file dropbox.com/s/gn9bzvkfbdi5j … .root?dl=0

6

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)

7

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)

8

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…

9

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");
  c->cd(0);
}

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?

10

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)

11

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

12

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).

13

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