void calc_eta_xsecs(){ double Eg[45]; double dEg[45]; double eff_eta_3piq[45]; double deff_eta_3piq[45]; double eff_eta_2g[45]; double deff_eta_2g[45]; double eff_eta_3pi0[45]; double deff_eta_3pi0[45]; double Nps[45]; double Neta_2g[45],dNeta_2g[45]; double Neta_3piq[45],dNeta_3piq[45]; double Neta_3pi0[45],dNeta_3pi0[45]; double regge[45]={0.445562,0.383408,0.335138,0.296947,0.266275,0.241265, 0.22051,0.202978,0.187936,0.174876,0.163436, 0.153342,0.14437,0.136341, 0.129108 ,0.122556, 0.116593,0.111144, 0.106146,0.101547,0.0973009, 0.09337,0.08972,0.08632,0.08314,0.08018,0.0774,0.07479,0.07234,0.07,0.06785, 0.06579,0.06384,0.06199,0.06024,0.05857,0.05699,0.05548,0.05405,0.05268, 0.05138,0.05013,0.04893,0.04779,0.04669}; //ifstream inps("ps_flux.dat"); ifstream inps("flux.dat"); for (int i=18;i<45;i++){ inps >> Nps[i]; } inps.close(); for (int i=0;i<45;i++){ Eg[i]=2.9+0.2*double(i); dEg[i]=0.1; } ifstream eff_2g_file("eta_2g_eff.dat"); for (int i=0;i<45;i++){ eff_2g_file >> eff_eta_2g[i]; eff_2g_file >> deff_eta_2g[i]; } eff_2g_file.close(); ifstream eff_3pi_file("eta_3piq_eff.dat"); for (int i=0;i<45;i++){ eff_3pi_file >> eff_eta_3piq[i]; eff_3pi_file >> deff_eta_3piq[i]; } eff_3pi_file.close(); ifstream eff_3pi0_file("eta_3pi_eff.dat"); for (int i=0;i<45;i++){ eff_3pi0_file >> eff_eta_3pi0[i]; eff_3pi0_file >> deff_eta_3pi0[i]; } eff_3pi0_file.close(); ifstream in("eta_2g_yields.dat"); for (int i=0;i<45;i++){ double dummy; in >> dummy; in >> Neta_2g[i]; in >> dNeta_2g[i]; in >> dummy; in >> dummy; } in.close(); ifstream in2("eta_3piq_yields.dat"); for (int i=0;i<45;i++){ double dummy; in2 >> dummy; in2 >> Neta_3piq[i]; in2 >> dNeta_3piq[i]; } in2.close(); ifstream in3("eta_6g_yields.dat"); for (int i=0;i<45;i++){ double dummy; in3 >> Neta_3pi0[i]; in3 >> dNeta_3pi0[i]; in3 >> dummy; in3 >> dummy; } in3.close(); double xsec_2g[45],xsec_3piq[45],xsec_3pi0[45]; double xsec_param[45],xsec_param_up[45],xsec_param_down[45]; double dxsec_2g[45],dxsec_3piq[45],dxsec_3pi0[45]; for (int i=0;i<45;i++){ double denom=Nps[i]; double eta_2g_numer=Neta_2g[i]; double r_eta_2g=eta_2g_numer/denom; double dr_eta_2g=0.; if (eta_2g_numer>0) dr_eta_2g=r_eta_2g*sqrt(1./denom+(dNeta_2g[i]*dNeta_2g[i]) /eta_2g_numer/eta_2g_numer+pow(0.2/39.41,2)); // printf("%f %f %f \n",1/denom,dNeta_2g[i],dNeta_2g_acc[i]); double scale_eta_2g=1e6/(1.293)/eff_eta_2g[i]/0.3941; //myflux << denom*1.249/(eff_ps[i]*1.65e-4) << endl; xsec_2g[i]=scale_eta_2g*r_eta_2g; dxsec_2g[i]=scale_eta_2g*dr_eta_2g; double eta_3piq_numer=Neta_3piq[i]; double r_eta_3piq=eta_3piq_numer/denom; double dr_eta_3piq=0.; if (eta_3piq_numer>0) dr_eta_3piq=r_eta_3piq*sqrt(1./denom+(dNeta_3piq[i]*dNeta_3piq[i]) /eta_3piq_numer/eta_3piq_numer +pow(0.28/22.92,2)); double scale_eta_3piq=1e6/(1.293)/eff_eta_3piq[i]/0.2292; xsec_3piq[i]=scale_eta_3piq*r_eta_3piq; dxsec_3piq[i]=scale_eta_3piq*dr_eta_3piq; double eta_3pi0_numer=Neta_3pi0[i]; double r_eta_3pi0=eta_3pi0_numer/denom; double dr_eta_3pi0=0.; if (eta_3pi0_numer>0) dr_eta_3pi0=r_eta_3pi0*sqrt(1./denom+(dNeta_3pi0[i]*dNeta_3pi0[i]) /eta_3pi0_numer/eta_3pi0_numer +pow(0.23/32.68,2)); double scale_eta_3pi0=1e6/(1.293)/eff_eta_3pi0[i]/0.3268; xsec_3pi0[i]=scale_eta_3pi0*r_eta_3pi0; dxsec_3pi0[i]=scale_eta_3pi0*dr_eta_3pi0; /* xsec_param[i]=0.288+28.2*pow(Eg[i],-2.11); double xsec_var=0.067*0.067+1.54*1.54*pow(Eg[i],-2.*2.11) +2.11*2.11*pow(Eg[i],2.*(-3.11))*0.4*0.4*28.2*28.2; //printf("%f\n",sqrt(xsec_var)); xsec_param_up[i]=xsec_param[i]+sqrt(xsec_var); xsec_param_down[i]=xsec_param[i]-sqrt(xsec_var); */ printf("%3.2f %f %f %f %f %f %3.2f+-%3.2f %3.2f+-%3.2f\n", Eg[i],Neta_2g[i],eff_eta_2g[i],Neta_3piq[i],eff_eta_3piq[i], denom,xsec_2g[i],dxsec_2g[i],xsec_3piq[i],dxsec_3piq[i] ); } TH2F *h=new TH2F("h","#eta cross section estimates",45,2.8,11.8,10,0,0.5); h->SetXTitle("E_{#gamma} [GeV]"); h->GetXaxis()->SetLabelSize(0.04); h->GetXaxis()->SetTitleSize(0.05); h->SetYTitle("#sigma [#mub]"); h->GetYaxis()->SetLabelSize(0.04); h->GetYaxis()->SetTitleSize(0.05); h->Draw(); TGraphErrors *g1=new TGraphErrors(45,Eg,xsec_2g,dEg,dxsec_2g); g1->SetMarkerStyle(20); g1->Draw("p"); TGraphErrors *g2=new TGraphErrors(45,Eg,xsec_3piq,dEg,dxsec_3piq); g2->SetMarkerColor(2); g2->SetLineColor(2); g2->SetMarkerStyle(20); g2->Draw("p"); TGraphErrors *g3=new TGraphErrors(45,Eg,xsec_3pi0,dEg,dxsec_3pi0); g3->SetMarkerColor(4); g3->SetLineColor(4); g3->SetMarkerStyle(20); g3->Draw("p"); TGraph *g5=new TGraph(45,Eg,regge); g5->Draw("c"); /* TGraph *g4=new TGraph(45,Eg,xsec_param_up); g4->SetLineStyle(2); g4->Draw("l"); TGraph *g5=new TGraph(45,Eg,xsec_param_down); g5->SetLineStyle(2); g5->Draw("l"); */ double oldxsec1[1]={0.462}; double doldxsec1[1]={0.128}; double Eold1[1]={4.}; double dEold1[1]={0.2}; double oldxsec2[1]={0.2}; double doldxsec2[1]={0.1}; double Eold2[1]={5.15}; double dEold2[1]={1.15}; TGraphErrors *g6=new TGraphErrors(1,Eold1,oldxsec1,dEold1,doldxsec1); g6->SetMarkerStyle(31); g6->Draw("p"); TGraphErrors *g7=new TGraphErrors(1,Eold2,oldxsec2,dEold2,doldxsec2); g7->SetMarkerStyle(30); g7->Draw("p"); TLegend *legend=new TLegend(0.45,0.6,0.85,0.85); legend->AddEntry(g1,"#eta#rightarrow2#gamma","p"); legend->AddEntry(g2,"#eta#rightarrow#pi^{+}#pi^{-}#pi^{0}","p"); legend->AddEntry(g3,"#eta#rightarrow3#pi^{0}","p"); legend->AddEntry(g5,"Laget Regge model","l"); legend->AddEntry(g6,"Bellenger,PRL21,1205","p"); legend->AddEntry(g7,"Struczinksi,NPB108,45","p"); legend->Draw(); }