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