--1556251094-491137576-904220092=:7405
Content-Type: TEXT/PLAIN; charset=US-ASCII
Hello Rooters,
I have some problems with histograms statistics after fitting.
I'm using ROOT 2.00/10 RedHat 5.0 with glibc.
At first below program (a part of bigger program)
fits some spectrum (sorry, a few minutes), write the results
(10 parameters of fitting as a histogram to the file), but
a few last errors of these parameters (ey[7],ey[8],ey[9])
are drawn uncorectly as I insert line
hout[p][r]->Fill(x[i],y[i]);
All is OK if I use the graph or
hout[p][r]->SetBinContent(i+1,y[i]);
I don't understand such behaviour. I prefer to use Fill,
because SetBinContent doesn't calculate Mean and RMS.
In both cases I have to calculate numbers of entries manually
and set it by TH1::SetEntries. I think there whould be nice to have
TH1::SetMean() and TH1::SetRMS() as well.
Next the graph hasn't got such statistics box. :-(
Nevertheless, I'd like to know whether is implemented for the
histogram in ROOT the calculation of the right mean and the RMS in case
such external errors?
Is it made by weights w[i]=1/ey[i]? Is it made automaticly?
I'm not sure if I insert it properly.
My fitting gives big errors as I use the Minos technique, so
I'm not sure how reliable my results are.
Nevertheless, I'm surprised that if I insert N_mean manually ( so
little diferent from the read value)
the results would change so dramaticly. Is my fitting so flexible?
Any comments are welcome.
Tadeusz
// spect.C
Int_t p_min=5;
Int_t p_max=5;
Int_t r_min=6;
Int_t r_max=6;
Int_t m_min=4;
Int_t m_max=32;
Int_t npar=10;
Int_t Nshow_max=500;
const Int_t M=35;
Double_t fac[M];
Double_t shift=-0.2;
Double_t step=0.05;
Double_t N_mean;
Int_t p;
Int_t r;
Int_t m;
Double_t Factorial(Int_t f)
{
Int_t b;
Double_t x;
if(f==0)
{
x=1.0;
return x;
}
x=1;
b=0;
do
{
b=b+1;
x=x*b;
}while(b!=f);
return x;
}
Double_t fitf(Double_t *x, Double_t *par)
{
Int_t N;
Double_t dN,lN;
Int_t k=x[0];
if(k>m) return 0;
Double_t kk=fac[k];
Double_t mm=fac[m];
Double_t mk=fac[m-k];
Double_t Binomial=mm/kk/mk;
Double_t a1,a2,a3,fun1;
Double_t At,Ad,a0;
At=TMath::Pi()*100*((r+1)*(r+1)-r*r);
Ad=3.24;
a0=Ad/At;
Double_t sum=0;
for(N=0;N<npar;N++)
{
lN=shift+N*step;
dN=N_mean*TMath::Power(10,lN);
a1=1-TMath::Exp(-dN*a0);
a2=TMath::Power(a1,k);
a3=TMath::Exp(-dN*(m-k)*a0);
fun1=Binomial*a2*a3;
sum=sum+par[N]*fun1;
}
return sum;
}
void spect()
{
TFile *fin=new TFile("dep.root");
TFile *fout = gROOT->FindObject("spect.root");
if(fout) fout->Close();
fout = new TFile("spect.root","RECREATE","ROOT file");
fin->cd();
TCanvas *page = new TCanvas("page");
page->Divide(1,2);
page->cd(1);
page_1->SetGrid();
page_1->SetLogy();
const Int_t P=35;
const Int_t R=20;
TH1S *hin[P][R];
TH1F *hout[P][R];
TF1 *fun[P][R];
for(Int_t i=0;i<M;i++)
fac[i]=Factorial(i);
char hname[20],hname1[20],hname2[20];
char dirname[50];
Int_t j,hit,np,sfit,Nshow,sum_Nshow,max_Nshow,mx;
Float_t x[npar],ex[npar],y[npar],ey[npar];
TF1 *func=new TF1("func",fitf,0,35,npar);
func->SetLineColor(2);
fout->cd();
TDirectory *din=fout->mkdir("din");
TDirectory *dout=fout->mkdir("dout");
for(p=p_min;p<=p_max;p++) {
for(r=r_min;r<=r_max;r++) {
sprintf(hname,"h%d",0);
TH2S *h0=(TH2S*)fin.Get(hname);
sprintf(hname1,"h%d",1);
TH1F *h1=(TH1F*)fin.Get(hname1);
sprintf(hname2,"h%d",2);
TH1F *h2=(TH1F*)fin.Get(hname2);
N_mean=h2->GetMean();
// N_mean=57.985262;
printf("N_mean=%f\n",N_mean);
for(Int_t i=0;i<npar;i++) {
y[i]=0;
ey[i]=0;
}
char pname[20];
sum_Nshow=0;
max_Nshow=0;
for(m=m_min;m<=m_max;m++) {
Nshow=h1->GetBinContent(m+1);
printf("oooo m=%d Nshow=%d oooo\n",m,Nshow);
if(Nshow>max_Nshow) {
max_Nshow=Nshow;
mx=m;
}
}
m=mx;
printf("********************************************************\n");
printf("p=%d r=%d m=%d max_Nshow=%d\n",p,r,m,max_Nshow);
printf("********************************************************\n");
if(max_Nshow>Nshow_max) {
din->cd();
sprintf(pname,"hin_p%d_r%d",p,r);
hin[p][r]=new TH1S(pname,"hin",m+2,0,m+1);
hin[p][r].SetMinimum(1e-2);
np=0;
for(j=0;j<=m;j++) {
hit=h0->GetCellContent(j+1,m+1);
hin[p][r]->Fill(j,hit);
np=np+hit;
}
hin[p][r].SetEntries(np);
for (Int_t ipar=0;ipar<npar;ipar++) {
func->SetParLimits(ipar,0,1e10);
func->SetParameter(ipar,1);
}
hin[p][r]->Fit("func","brew","",0,m);
fun[p][r]=hin[p][r]->GetFunction("func");
for(Int_t i=0;i<npar;i++) {
y[i]=fun[p][r]->GetParameter(i);
ey[i]=fun[p][r]->GetParError(i);
}
sum_Nshow=sum_Nshow+Nshow;
} else
for(Int_t i=0;i<npar;i++) {
y[i]=0;
ey[i]=0;
}
page->cd(2);
page_2->SetGrid();
page_2->SetLogy();
dout->cd();
hout[p][r]=new TH1F("hout","hout",npar,shift-step,shift+(npar+1)*step);
hout[p][r].SetLineWidth(2);
hout[p][r].SetLineStyle(1);
hout[p][r].SetLineColor(1);
hout[p][r].SetXTitle("dN");
hout[p][r].SetYTitle("N");
hout[p][r].SetMinimum(1e-5);
hout[p][r].SetMaximum(1e4);
Int_t nm;
Float_t Nmi,dNmi;
sfit=0;
// hout[p][r]->Sumw2();
for(Int_t i=0;i<npar;i++)
{
x[i]=shift+i*step;
ex[i]=0.001;
hout[p][r]->Fill(x[i],y[i]);
// hout[p][r]->SetBinContent(i+1,y[i]);
hout[p][r]->SetBinError(i+1,ey[i]);
printf("y[%d]=%f ey[%d]=%f\n",i,y[i],i,ey[i]);
sfit=sfit+y[i];
}
sfit=(int) sfit;
hout[p][r]->SetEntries(sfit);
hout[p][r]->Draw("A");
char hn[20];
sprintf(hn,"hout_p%d_r%d",p,r);
hout[p][r].SetName(hn);
}
}
fout.Write();
delete fin;
delete fout;
}
-- Tadeusz Pytlos mailto:pytlos@fizwe5.fic.uni.lodz.pl Lodz, Poland--1556251094-491137576-904220092=:7405 Content-Type: APPLICATION/octet-stream; name="dep.root" Content-Transfer-Encoding: BASE64 Content-ID: <Pine.LNX.3.96.980827141451.7405B@fizwe5.fic.uni.lodz.pl> Content-Description:
cm9vdAAATioAAABAAAAGrAAAAJQAAAA/AAAAAgAAADYEAAAAAQAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAFQAAgAAACkONr6lACsAAQAAAEAAAAAA BVRGaWxlCWRlcG4ucm9vdAAJZGVwbi5yb290AAABDja+pQ42wG4AAAC6AAAA NgAAAEAAAAAAAAAF8gAAAD8AAgAAABQONsBvACsAAQAAAJQAAABABVRGaWxl CWRlcG4ucm9vdAAAAQAAANMAAAD3AAEAAAasdzWUAP///9sAwwAAAEAFVEZp bGUJZGVwbi5yb290AAABAAAA+Hc1lAAAAAJmAAIAAAxIDjbAVAAzAAEAAAD4 AAAAQARUSDJTAmgwEGstbSBkaXN0cmlidXRpb25DUwgqAgBIDADdlj9oU1EU xr97X14SkhdTtEHQDrUVKqJSwVZo1DRCtw62UKE4tKYdTItVTJWqSzdBBQVx dhDN4tCpc1C3QEsXEyelg9hNoV3Fc+7NX8yQvJPJ3+Pm3Pfyzsc55513EihU D8KhpW8N9yyfvd2/mM2t3stm7q9m76zU7qG7nD27S62zhztLTua69XbXbq5l c6QCqD/miyvJt6XI5cFS5NLBxxfsxHu6c5AWrnr82UCz2MNuij1qV8y42PQa cX6Of9sq8I4sm7otfrG2vLPO9vPT72Y1Y2RNxf53uO7d4TSCdiNFpTBi2kCM 8wxFHIe2ZxJCT9QuLiKAWm/4xRtzctjECdISKsV+BzL4gWNGSaQVDutXeIcB hKS1CuX1e+QxREquLCZ3W4/jNc5IdYDANbWJuxQTPz0RTlRNYB5HK9nZwxf6 q7qAaZwkJVbQfnXI9Y0awAL6KCLhs9OzSmEFR0iJO6pZraP49DmVxBwO16rU oX8d3asO4QZV3DE1EiipDexjipQiZrY0Z9cR6jpeIo2EeXY2Jr9R7eE5RhGv ZOdXhfmEJZopMfkswGPM4BQ8ikhQI8MDnKcqxU1urWNqffVfxjBM3RSG/bVq 16sVaXpTvMpbV6Vx3z4j6KHceGL686+TQC+iYhUmRoedl9XFVG0nRCu9JO/M IEXkmndOikcdEKR45P9ZeMJxTlIdniNcpW7AKlZJpqdSv+wmnTj4YOz06FSB bXCyzOd/AQAAAT8AAgAAAgQONsBcADEAAQAAA14AAABABFRIMUYCaDEObSBk aXN0cmlidXRpb25DUwgFAQAEAgBjYGSAQCBgBmKmDEO+XIWUzOKSosyk0pLM /DyoPFgN80sIy74BpJ5BFagFLA7Ry1qRWJFZDDSDgYHxH1jCznrxdS5b5etc Nl8P9YM0gdhAlcpAzODEAyKRAKphlcQaBtYCdRECoBpWRZlhDMwvHO6fOwBi AWkQBacdRT9/ANNBRr4NQPpI1wMwRgVgYxlAIYYEHBqAgeDAwOC8u4HBpS+B wVUqgMG1xYPBddYJBlezDQyuRg0MrrwJDC5XFBhcZj1gcJl4gMGlA0hXMTC4 VADZ6UDaDUgLHmBwftDA4LwBiFuA5gUAMSfQ7CdAvAeIQxjAge0kAHTmDCDN AKQNgPY/AGKg/QAAAAFVAAIAAAH5DjbAYAAnAAEAAASdAAAAQARUSDFGAmgy BGhvdXRDUwglAQD5AQBjYGSAQCBgZmDgYMowYsnILy2BijKBZJhfQnj2DSBV DDxAhSA+A1gHA2tFYkVmMVOKH1AxA+M/sJSd9eLrXLbK17lsvh7qB2kDsYFq uQ4qMDC4VDGgAlTjKkHGMRJjGlgP1FEIgGpaFcg0BrINY2B+4bBZA8wC0nde b3FhdUxvPbU/05bLUTD6uYDixc2O1TY7DggdTnI/HMBg1SyUD1YMATwQCqh3 y+KLrLU3GRy2X11jfU+SwcGJ9VPuqv0KDvKXb/lpVDI4CHqwq2+dwODA+3+3 /Bk7BQehqftZi6IcHBRev797N1fBwS5FxHvJTohhKADsbAaoVa4MD0tcDVdu MZq5WtC8pHyb2Y67GiaLt7YY3G8pM7ffFm/Gdr3STnPpJ5BaAAAAALoAAgAA AI8ONsBuACsAAQAABfIAAABABVRGaWxlCWRlcG4ucm9vdAAAAAADAAACZgAC AAAMSA42wFQAMwABAAAA+AAAAEAEVEgyUwJoMBBrLW0gZGlzdHJpYnV0aW9u AAABPwACAAACBA42wFwAMQABAAADXgAAAEAEVEgxRgJoMQ5tIGRpc3RyaWJ1 dGlvbgAAAVUAAgAAAfkONsBgACcAAQAABJ0AAABABFRIMUYCaDIEaG91dA== --1556251094-491137576-904220092=:7405--