Stat & Fit box handling

From: Tadeusz Pytlos (pytlos@fizwe5.fic.uni.lodz.pl)
Date: Wed Oct 21 1998 - 14:59:36 MEST


Dear Rooters,
I have problems to receive full information
about my histos, I mean stat box and fit box together.
The example shows that I succesed only with one.
How receive stat&fit box in all graphes?
Thank you in advance,
                 Tadeusz
P.S. Sorry for the long example.

--
Tadeusz Pytlos        
mailto:pytlos@fizwe5.fic.uni.lodz.pl 
Lodz, Poland                                                 

\\ gama1.C

Int_t p=15;
Int_t r=6;
Int_t width=100;
Int_t shift=width/2;
Int_t npar=10;
Int_t MM=35;
Int_t m=12;
Double_t Ns1=3000;
Double_t Ns2,Ns3,Ns4;
Double_t fac[MM+1];
// 1
// Double_t pp[]={0,1,0};
// 2
// Double_t pp[]={0.5,0,0.5};
// 3
Double_t pp[]={0.2,0.4,0.4};
Double_t N_mean[]={100,400,630};
Double_t M[]={5,8,100};

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 Gamma(Double_t z)
{
  Double_t x,f,g,gamma;
  Double_t p2pil=0.91893853;
  if(z<=0) exit;
  x=z;
  f=1.0;
  if(z>6.0)  goto L100;
  x=6.0+z;
  f=1.0/z;
  for(Int_t i=1;i<6;i++) f=f/(z+float(i));
L100:
  g=p2pil+(x-0.5)*TMath::Log(x)-x;
  g=g+TMath::Log(1.0+1.0/12.0/x+1.0/288.0/TMath::Power(x,2)
    -139.0/51840.0/TMath::Power(x,3));
  gamma=TMath::Exp(g)*f;
  return gamma;
}

Double_t fitf(Double_t *x, Double_t *par)
{
  Int_t N,dN;
  Int_t k=x[0]+0.5;
  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;
  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++)
  {
    dN=shift+width*N;
    a1=1-TMath::Exp(-dN*a0);
    a2=TMath::Power(a1,k);
    a3=TMath::Exp(-dN*(m-k)*a0);
    sum=sum+par[N]*Binomial*a2*a3;
  }
  return sum; 
}

Double_t fGamma(Double_t *x,Double_t *par)
{
  Double_t fun1,fun2,fun3;
  Double_t Nm=x[0];
  Double_t b1=TMath::Power(par[4]/par[3],par[4]);
  Double_t b2=TMath::Power(Nm,par[4]-1);
  Double_t b3=TMath::Exp(-par[4]*Nm/par[3]);
  fun1=b1*b2*b3/Gamma(par[4]);
  if((par[5]>0) && (par[6]>0)) {
  Double_t d1=TMath::Power(par[6]/par[5],par[6]);
  Double_t d2=TMath::Power(Nm,par[6]-1);
  Double_t d3=TMath::Exp(-par[6]*Nm/par[5]);
  fun2=d1*d2*d3/Gamma(par[6]);
  } else fun2=0;
  if((par[7]>0) && (par[8]>0)) {
  Double_t e1=TMath::Power(par[8]/par[7],par[8]);
  Double_t e2=TMath::Power(Nm,par[8]-1);
  Double_t e3=TMath::Exp(-par[8]*Nm/par[7]);
  fun3=e1*e2*e3/Gamma(par[8]);
  } else fun3=0;
  Double_t fun123=Ns4*width*(par[0]*fun1+par[1]*fun2+par[2]*fun3);
  return fun123;
}

Double_t hGamma(Int_t k,Double_t N_mean,Double_t M)
{
  Double_t kk=fac[k];
  Double_t mm=fac[m];
  Double_t mk=fac[m-k];
  Double_t Bi1=mm/kk/mk;
  Double_t At,Ad,a0,sumc,sumd;
  At=TMath::Pi()*100*((r+1)*(r+1)-r*r);
  Ad=3.24;
  a0=Ad/At;
  Double_t c1=TMath::Power(M/N_mean,M);
  sumc=0;
  for(Int_t i=0;i<=k;i++) {
    Double_t ii=fac[i];
    Double_t ik=fac[k-i];
    Double_t Bi2=kk/ii/ik;
    Double_t zn=TMath::Power(-1,i);
    Double_t mian=(m-k)*a0+M/N_mean+i*a0;
    Double_t c2=TMath::Power(mian,-M);
    sumc=sumc+Bi2*zn*c2;
  }
  sumd=Bi1*c1*sumc;
  return sumd;
}

void gama1() 
{
   page = new TCanvas("page","Probability");
   page->SetFillColor(10);
   page->Divide(2,2);
   gStyle->SetOptStat(1111);
   gStyle->SetOptFit(111);


   Int_t kpad,ans;
   char padname[20];

// pad1-------------------------------------------------
   kpad=0;
   sprintf(padname,"page_%d",kpad+1);
   TPad *pad = (TPad*)page->GetPrimitive(padname);
   pad->cd();
   pad->SetGrid();
   pad->GetFrame()->SetFillColor(42);
   pad->GetFrame()->SetBorderMode(-1);
   pad->GetFrame()->SetBorderSize(5);
   pad->Draw();

   for(Int_t i=0;i<=MM;i++)
     fac[i]=Factorial(i);

   TH1F *h1=new TH1F("h1","h1",npar,shift-width/2,
                       shift-width/2+npar*width);
   h1.SetLineWidth(2);
   h1.SetXTitle("dN");
   h1.SetYTitle("N");
   h1.SetMinimum(0);

   TF1 *f1=new TF1("f1",fGamma,shift-width/2,shift-width/2+npar*width,9);
   f1->SetLineStyle(1);
   f1->SetLineColor(4);
   f1->SetParameter(0,pp[0]);
   f1->SetParameter(1,pp[1]);
   f1->SetParameter(2,pp[2]);
   f1->SetParameter(3,N_mean[0]);
   f1->SetParameter(4,M[0]);
   f1->SetParameter(5,N_mean[1]);
   f1->SetParameter(6,M[1]);
   f1->SetParameter(7,N_mean[2]);
   f1->SetParameter(8,M[2]);
   f1->Draw();

   page.Update();
// pad2---------------------------------------------------------
   sprintf(padname,"page_%d",kpad+2);
   TPad *pad = (TPad*)page->GetPrimitive(padname);
   pad->cd();
   pad->SetGrid();
   pad->GetFrame()->SetFillColor(42);
   pad->GetFrame()->SetBorderMode(-1);
   pad->GetFrame()->SetBorderSize(5);
   pad->Draw();

   TH1F *h2=new TH1F("h2","h2",m+1,-0.5,m+0.5);   
   h2.SetLineWidth(2);
   h2.SetLineStyle(1);
   h2.SetLineColor(1);
   h2.SetXTitle("k");
   h2.SetYTitle("N");
   h2.SetMinimum(0);

   Float_t w2;

   Ns2=0;
   for(Int_t k=0;k<=m;k++) {
      w2=Ns1*(pp[0]*hGamma(k,N_mean[0],M[0])+pp[1]*hGamma(k,N_mean[1],M[1])
                 +pp[2]*hGamma(k,N_mean[2],M[2]));
      h2->Fill(k,w2);
      Ns2=Ns2+w2;
   }
   printf("Ns2=%f\n",Ns2);
   Ns2=(int) Ns2;
   h2->SetEntries(Ns2);
   h2.Draw();


   TH1F *h22=new TH1F("h22","h22",m+1,-0.5,m+0.5);   
//   gRandom->SetSeed(10);
   h22->FillRandom(h2,Ns2);
   h22.SetLineColor(2);
   h22->Draw("same");

   page.Update();
// pad3----------------------------------------------------------
   sprintf(padname,"page_%d",kpad+3);
   TPad *pad = (TPad*)page->GetPrimitive(padname);
   pad->cd();
   pad->SetGrid();
   pad->GetFrame()->SetFillColor(42);
   pad->GetFrame()->SetBorderMode(-1);
   pad->GetFrame()->SetBorderSize(5);
   pad->Draw();

   TH1F *h3=h22->Clone();
   h3.SetLineWidth(2);
   h3.SetLineStyle(1);
   h3.SetLineColor(1);
   h3.SetXTitle("k");
   h3.SetYTitle("N");
   h3.SetMinimum(0);
   
   TF1 *f2=new TF1("f2",fitf,0,m,npar);
   f2->SetLineStyle(1);
   f2->SetLineColor(7);
   for (Int_t ipar=0;ipar<npar;ipar++) {
      f2->SetParLimits(ipar,0,1.2*Ns1);
      f2->SetParameter(ipar,0.1*Ns1);
   }

   new TMinuit(npar);
   h3->Fit("f2","b");
   TF1 *ff2=gROOT->GetFunction("f2");
// pad4----------------------------------------------------------
   sprintf(padname,"page_%d",kpad+4);
   TPad *pad = (TPad*)page->GetPrimitive(padname);
   pad->cd();
   pad->SetGrid();
   pad->GetFrame()->SetFillColor(42);
   pad->GetFrame()->SetBorderMode(-1);
   pad->GetFrame()->SetBorderSize(5);
   pad->Draw();

   TH1F *h4=new TH1F("h4","h4",npar,shift-width/2,
                       shift-width/2+npar*width);
   h4.SetLineWidth(2);
   h4.SetXTitle("dN");
   h4.SetYTitle("N");
   h4.SetMinimum(0);
   h4.SetMaximum(0.5*Ns1);

   Float_t x4[npar],y4[npar],ex4[npar],ey4[npar];

   Ns4=0;
   for(Int_t i=0;i<npar;i++)
   {
     x4[i]=shift+i*width;
     y4[i]=ff2->GetParameter(i);
     ex4[i]=width/2;
     ey4[i]=ff2->GetParError(i);
     h4->Fill(x4[i],y4[i]);
     Ns4=Ns4+y4[i];
   }
   Ns4=(int) Ns4;
   h4->SetEntries(Ns4);
   h4->Draw("A");

   TF1 *fgam=new TF1("fgam",fGamma,shift-width/2,
                           shift-width/2+npar*width,9);
   fgam->SetLineColor(2);
   for (Int_t ipar=0;ipar<3;ipar++) {
      fgam->SetParLimits(ipar,0,1);
   }
   fgam->SetParLimits(3,1e-3,npar*width);
   fgam->SetParLimits(4,1e-3,20);
   fgam->SetParLimits(5,1e-3,npar*width);
   fgam->SetParLimits(6,1e-3,20);
   fgam->SetParLimits(7,1e-3,npar*width);
   fgam->SetParLimits(8,1e-3,20);

   fgam->SetParameter(0,pp[0]);
   fgam->SetParameter(1,pp[1]);
   fgam->SetParameter(2,pp[2]);
   fgam->SetParameter(3,N_mean[0]);
   fgam->SetParameter(4,M[0]);
   fgam->SetParameter(5,N_mean[1]);
   fgam->SetParameter(6,M[1]);
   fgam->SetParameter(7,N_mean[2]);
   fgam->SetParameter(8,M[2]);

   TGraphErrors *gr4=new TGraphErrors(npar,x4,y4,ex4,ey4);
   gr4->SetTitle("TGraphErrors Example");
   gr4->SetMarkerStyle(21);
   gr4->SetMarkerSize(0.7);
   gr4->SetMarkerColor(3);
   gr4->Draw("P");
   gr4->Fit("fgam","br");
   f1->Draw("same");

   page.Update();

   TPostScript mps("gama.eps",113);
   page.Draw();
   mps.Close();
}


--
Tadeusz Pytlos        
mailto:pytlos@fizwe5.fic.uni.lodz.pl 
Lodz, Poland                                                 



This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:34:38 MET