Tadeusz,
You have a mistake in your line creating the hout histogram
As a consequence you fill the same bin twice and your errors
are messed.
Change the line:
hout[p][r]=new
TH1F("hout","hout",npar,shift-step,shift+(npar+1)*step);
to
hout[p][r]=new TH1F("hout","hout",npar,shift-step,shift+npar*step);
Rene Brun
On Thu, 27 Aug 1998, Tadeusz Pytlos wrote:
> 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
>
This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:34:37 MET