analyse_phi.C File Reference

Go to the source code of this file.

Functions
void	confgraph (TGraph *g, TString tit, int col=1, int marker=20)
TH1F *	createHistoGraph (TGraph *g, TString tit="", double xmin=0, double xmax=0, double ymaxf=1.05)
void	analyse_phi (double t=0.5, double sigs=100, int steps=25, double emin=2.25, double emax=2.35, double m=2.3, double G=0.015)

Function Documentation

void analyse_phi	(	double	t = 0.5,
		double	sigs = 100,
		int	steps = 25,
		double	emin = 2.25,
		double	emax = 2.35,
		double	m = 2.3,
		double	G = 0.015
	)

Definition at line 50 of file analyse_phi.C.

References c1, confgraph(), createHistoGraph(), dx, dy, effs, f1, i, m, printf(), sqrt(), t, x, and y.

{
        //gStyle->SetStatX(0.5);
        //gStyle->SetStatY(0.40);
        gStyle->SetStatW(0.20);
        gStyle->SetStatH(0.12);
        
        // C. Evangelista et al., Phys. Rev. D57, 5370 (1998).
        double x[25]  = {2.219, 2.221, 2.221, 2.222, 2.224, 2.226, 2.226, 2.228, 2.229, 2.231, 2.231, 2.233, 2.235, 2.236, 2.242, 2.247, 2.254,  2.255, 2.272, 2.289, 2.307, 2.333, 2.360, 2.404, 2.430};
        double dx[25] = {0.0};
        
        double y[25]  = {4.15,  4.16,  3.87,  3.94,  4.40,  3.90,  3.68,  3.00,  3.20,  3.73,  2.97,  4.20,  3.52,  2.48,  3.38,  2.68,  2.63,   3.01,  2.25,  2.15,  1.97,  2.04,  2.04,  1.61,  1.93 };
        double dy[25] = {0.34,  0.36,  0.54,  0.35,  0.36,  0.34,  0.34,  0.30,  0.32,  0.35,  0.33,  0.35,  0.34,  0.30,  0.45,  0.39,  0.35,   0.16,  0.17,  0.21,  0.24,  0.14,  0.12,  0.18,  0.11 };


        TGraphErrors *gphi = new TGraphErrors(25,x,y,dx,dy);
        
        confgraph(gphi,"#phi#phi cross section [JetSet];E_{cm} [GeV];#sigma [#mub]");
        TH1F *hphi = createHistoGraph(gphi,"#phi#phi cross section [JetSet];E_{cm} [GeV];#sigma [#mub]",2.21,2.435);
        
        //TF1 *f1=new TF1("f1","expo",2.21,2.435);
        TF1 *f1=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",2.21,2.435);
        f1->SetParameters(3,2.23,0.01,1);       

        // inelasic cross section in mb
        TF1 *finel=new TF1("finel","[0]+[1]/x^[2]-[3]/x^[4]");
        finel->SetParameters(32,82,0.67,45,0.85);
        finel->SetRange(1.,100.);
        
        TCanvas *c1=new TCanvas("c1","c1",20,20,1500,500);
        c1->Divide(3,1,0.0001,0.0001);
        
        c1->cd(1);
        hphi->Draw();
        gphi->Fit("f1","R");
        gphi->Draw("Psame");

        TLine l;
        l.SetLineColor(4);
        l.SetLineStyle(9);
        l.SetLineWidth(2);
        l.DrawLine(emin,0,emin,hphi->GetMaximum()*0.5); 
        l.DrawLine(emax,0,emax,hphi->GetMaximum()*0.5); 
        
        
        double effs   = 20322./100000.;     // signal efficiency: cut = abs(summ4-2.15)<0.15&&abs(f4cxd0m-1.02)<0.012&&abs(f4cxd1m-1.02)<0.012&chi24c<50        
        double effg   = 35./840.5e6;        // generic background efficiency
        //double sigs   = 100.;               // signal cross section [nb]
        double fBR    = 0.489*0.489;        // BR(phi->K+ K-)^2
        
        //double effs   = 22463./100000.;     // signal efficiency: cut = abs(summ4-2.15)<0.15&&abs(f4cxd0m-1.02)<0.012&&abs(f4cxd1m-1.02)<0.012&chi24c<50      
        //double effg   = 77./840.5e6;        // generic background efficiency
        
        double sigphi = f1->Eval(2.3)*1000; // phi phi cross section [nb]
        double sigg   = 60e6;               // generic cross section [nb]
        
        //double t      = 0.5;                // [d]
        double Lint   = 788;               // [1./(nb*d)]
        
        double S      = sigs*effs*fBR*t*Lint;    // expected number of signals
        double Bphi   = sigphi*effs*fBR*t*Lint;  // expected number of non-res phi phi       
        double Bg     = sigg*effg*t*Lint;        // expected number of 
        
        cout <<"S="<<S<<"  B_phi(sig ="<<sigphi<<" nb)="<<Bphi<<"  B_g="<<Bg<<endl;


        c1->cd(2);
        
        // signal cross section (BW)
        TF1 *fbw = new TF1("fbw","[0]*[3]*TMath::BreitWigner(x,[1],[2])");
        fbw->SetParameters(1,m,G,1);
        fbw->SetParNames("#sigma_{S}","m","#Gamma","f_{scale}");

        // background function
        TF1 *fbg=new TF1("fbg","pol2",2.21,2.435);

        // combined fit function
        TF1 *ffit = new TF1("ffit","[0]*TMath::BreitWigner(x,[1],[2]) + pol2(3)");
        ffit->SetParNames("A","m","#Gamma","a_{0}","a_{1}","a_{2}","a_{3}");
        
        double stsize = (emax-emin)/(steps-1);
        
        TGraphErrors *gscan = new TGraphErrors(steps);
        TGraphErrors *gscansig = new TGraphErrors(steps);
        TGraphErrors *gscanbg = new TGraphErrors(steps);
        confgraph(gscan,Form("#phi#phi scan (t=%.2fd);E_{cm} [GeV];counts",t));
        confgraph(gscanbg,Form("#phi#phi scan (t=%.2fd);E_{cm} [GeV];counts",t),16);
        confgraph(gscansig,Form("#phi#phi scan (t=%.2fd);E_{cm} [GeV];counts",t));
        
        double st_t = t/steps;
                
                
                
        for (int i=0;i<steps;++i)
        {
                double ecm = emin + i*stsize;
                
                double st_sigs   = fbw->Eval(ecm)/fbw->Eval(m)*sigs;
                double st_sigphi = f1->Eval(ecm)*1000; // mub -> nb
                double st_sigg   = finel->Eval(ecm)*1e6; // mub -> nb
                
                double yieldS    = gRandom->Poisson(st_sigs * effs * fBR * st_t * Lint);
                double yieldBphi = gRandom->Poisson(st_sigphi * effs * fBR * st_t * Lint);
                double yieldBg   = gRandom->Poisson(st_sigg * effg * st_t * Lint);
                
                printf("(%2d) Ecm = %5.3f GeV : sig_S = %.2f nb (S = %6.0f)   sig_phi = %.2f nb (B = %6.0f)    sig_gen = %.2f mub (Bg = %5.0f)\n", i, ecm,  st_sigs, yieldS, st_sigphi, yieldBphi, st_sigg/1000, yieldBg);

                double yield     = (yieldS + yieldBg + yieldBphi);
                
                gscan->SetPoint(i,ecm, yield);
                gscan->SetPointError(i, 0, sqrt(yield));
                
                gscanbg->SetPoint(i,ecm, gRandom->Poisson(yieldBg));
                gscanbg->SetPointError(i, 0, sqrt(yieldBg));
        }       
                
        ffit->SetParameters(sigs,m,G,1,1,1);
        ffit->FixParameter(1,m);
        ffit->FixParameter(2,G);
        
        gscan->Fit("ffit");
        ffit->ReleaseParameter(1);
        ffit->ReleaseParameter(2);
        gscan->Fit("ffit","m");
        
        fbg->SetParameters(ffit->GetParameter(3),ffit->GetParameter(4),ffit->GetParameter(5),ffit->GetParameter(6));

        TH1F *hscan = createHistoGraph(gscan,Form("#phi#phi scan (t = %.2fd, #sigma = %.0fnb);E_{cm} [GeV];counts",t,sigs),0,0,1.3);
        hscan->Draw();  
        gscan->Draw("Psame");   
        gscanbg->Draw("Psame");
        
        for (int i=0;i<steps;++i)
        {
                double x,y,dx,dy;
                gscan->GetPoint(i,x,y);
                gscansig->SetPoint(i,x,y-fbg->Eval(x));
                gscansig->SetPointError(i,gscan->GetErrorX(i),gscan->GetErrorY(i));
                
                //cout <<x<<" "<<fbg->Eval(x)<<endl;
        }
        
        c1->cd(3);
        TH1F *hscansig = createHistoGraph(gscansig,Form("#phi#phi scan (t = %.2fd, #sigma = %.0fnb);E_{cm} [GeV];counts",t,sigs),0,0,1.3);
        hscansig->SetMinimum(-100);
        hscansig->Draw();
        gscansig->Draw("Psame");
        fbw->FixParameter(3,1); 
        gscansig->Fit("fbw");
        fbw->SetParameter(0,1);
        fbw->FixParameter(3, 1./fbw->GetMaximum(emin,emax)*(effs*fBR*st_t*Lint) );
        fbw->SetParameter(0,TMath::MaxElement(gscansig->GetN(),gscansig->GetY())/(effs*fBR*st_t*Lint));
        gscansig->Fit("fbw");
        
        
        c1->Update();
        TPaveStats *st = (TPaveStats *)gscan->GetListOfFunctions()->FindObject("stats");
        if (st)
        {
                cout <<"good"<<endl;
                st->SetX1NDC(0.5); //new x start position
                st->SetX2NDC(0.5); //new x end position
                c1->Modified();
                c1->Update();
        }
        else cout <<"***"<<endl;
                
}

void confgraph	(	TGraph *	g,
		TString	tit,
		int	col = 1,
		int	marker = 20
	)

Definition at line 3 of file analyse_phi.C.

References col.

{
        g->GetHistogram()->SetTitle(tit);
        g->GetHistogram()->SetMinimum(0);
        g->SetLineColor(col);
        g->SetMarkerColor(col);
        g->SetMarkerStyle(marker);
        g->SetMarkerSize(1);
        g->SetLineWidth(2);
}

TH1F* createHistoGraph	(	TGraph *	g,
		TString	tit = "",
		double	xmin = 0,
		double	xmax = 0,
		double	ymaxf = 1.05
	)

Definition at line 16 of file analyse_phi.C.

References cnt, h, xmax, and xmin.

{
        static int cnt=0;
        if (tit=="") tit=g->GetTitle();
        
        double ymax  = TMath::MaxElement(g->GetN(), g->GetY());
        double dymax = g->GetErrorY(TMath::LocMax(g->GetN(), g->GetEY()));
        
        if (xmin==xmax)
        {
                xmin = TMath::MinElement(g->GetN(), g->GetX());
                xmax = TMath::MaxElement(g->GetN(), g->GetX());
        }       
        TH1F *h=new TH1F(Form("h%03d",cnt++),tit,g->GetN(),xmin-0.5*(xmax-xmin)/g->GetN(), xmax+0.5*(xmax-xmin)/g->GetN());
        
        h->GetXaxis()->SetLabelSize(0.05);
        h->GetXaxis()->SetTitleSize(0.05);
        h->GetYaxis()->SetLabelSize(0.05);
        h->GetYaxis()->SetTitleSize(0.05);
        h->GetYaxis()->SetTitleOffset(1.4);
        
        h->GetXaxis()->SetNdivisions(505);
        
        h->SetMaximum((ymax+dymax)*ymaxf);
        h->SetMinimum(0);
        
        h->SetStats(0);
                
        return h;
}