void calc_phi_xsecs(){
  double Eg[45];
  double dEg[45];

  double eff_phi_3pi[45]; 
  double deff_phi_3pi[45];
  double eff_phi_3g[45];
  double deff_phi_3g[45];
  double eff_phi_kpkm[45];
  double deff_phi_kpkm[45];

  double Nps[45];
  double Nphi_3g[45],dNphi_3g[45];
  double Nphi_3pi[45],dNphi_3pi[45];
  double Nphi_kpkm[45],dNphi_kpkm[45];

  //ifstream inps("ps_flux.dat"); 
  ifstream inps("flux.dat");
  for (int i=18;i<45;i++){
    inps >> Nps[i];
    //Nps[i]=2e11;
  }
  inps.close();

  for (int i=0;i<45;i++){
    Eg[i]=2.9+0.2*double(i);
    dEg[i]=0.1;
  }
  ifstream eff_3g_file("phi_3g_eff.dat");
  for (int i=0;i<45;i++){
    eff_3g_file >> eff_phi_3g[i];
    eff_3g_file >> deff_phi_3g[i];
  }
  eff_3g_file.close(); 

  ifstream eff_3pi_file("phi_3pi_eff.dat");
  for (int i=0;i<45;i++){
    eff_3pi_file >> eff_phi_3pi[i]; 
    eff_3pi_file >> deff_phi_3pi[i];
  }
  eff_3pi_file.close();

  ifstream eff_kpkm_file("phi_kpkm_eff.dat");
  for (int i=0;i<45;i++){
    eff_kpkm_file >> eff_phi_kpkm[i]; 
    eff_kpkm_file >> deff_phi_kpkm[i];
    //    eff_phi_kpkm[i]=0.1;
  
  }
  eff_kpkm_file.close();

  ifstream in2("phi_kpkm_yields.dat");
  for (int i=0;i<45;i++){
    double dummy;
    in2 >> dummy;
    in2 >> Nphi_kpkm[i];
    in2 >> dNphi_kpkm[i];
  }
  in2.close(); 

  ifstream in("phi_yields.dat");
  for (int i=0;i<45;i++){
    double dummy;
    in >> dummy;
    in >> Nphi_3g[i];
    in >> dNphi_3g[i];
    in >> Nphi_3pi[i];
    in >> dNphi_3pi[i];
  }
  in.close(); 

  double xsec_3g[45],xsec_3pi[45],xsec_kpkm[45];
  double xsec_param[45],xsec_param_up[45],xsec_param_down[45];
  double dxsec_3g[45],dxsec_3pi[45],dxsec_kpkm[45];

  for (int i=0;i<45;i++){
    double denom=Nps[i];
    double phi_3g_numer=Nphi_3g[i];
    double r_phi_3g=phi_3g_numer/denom;
    double dr_phi_3g=0.;
    if (phi_3g_numer>0)
      dr_phi_3g=r_phi_3g*sqrt(1./denom+(dNphi_3g[i]*dNphi_3g[i])
		       /phi_3g_numer/phi_3g_numer);
    // printf("%f %f %f \n",1/denom,dNphi_3g[i],dNphi_3g_acc[i]);
    
    double scale_phi_3g=1e6/(1.293)/eff_phi_3g[i]
      /0.01304;
    
    xsec_3g[i]=scale_phi_3g*r_phi_3g;
    dxsec_3g[i]=scale_phi_3g*dr_phi_3g;

    double phi_3pi_numer=Nphi_3pi[i];
    double r_phi_3pi=phi_3pi_numer/denom;
    double dr_phi_3pi=0.;
    if (phi_3pi_numer>0)
      dr_phi_3pi=r_phi_3pi*sqrt(1./denom+(dNphi_3pi[i]*dNphi_3pi[i])
		       /phi_3pi_numer/phi_3pi_numer);
    double scale_phi_3pi=1e6/(1.293)/eff_phi_3pi[i]
      /0.1525;

    xsec_3pi[i]=scale_phi_3pi*r_phi_3pi;
    dxsec_3pi[i]=scale_phi_3pi*dr_phi_3pi;  

    double phi_kpkm_numer=Nphi_kpkm[i];
    double r_phi_kpkm=phi_kpkm_numer/denom;
    double dr_phi_kpkm=0.;
    if (phi_kpkm_numer>0)
      dr_phi_kpkm=r_phi_kpkm*sqrt(1./denom+(dNphi_kpkm[i]*dNphi_kpkm[i])
		       /phi_kpkm_numer/phi_kpkm_numer);
    double scale_phi_kpkm=1e6/(1.293)/eff_phi_kpkm[i]
      /0.492;
  cout << eff_phi_kpkm[i] << endl;

    xsec_kpkm[i]=scale_phi_kpkm*r_phi_kpkm;
    dxsec_kpkm[i]=scale_phi_kpkm*dr_phi_kpkm;

    xsec_param[i]=0.378+0.00117*pow(Eg[i],0.94)+0.0404*log(Eg[i]);
    double xsec_var=0.0309*0.0309+0.007*0.007*pow(Eg[i],2.*0.94)
      +0.00117*0.00117*0.94*0.94*pow(Eg[i],2.*(-0.06))*3.47*3.47
      +0.017*0.017*log(Eg[i])*log(Eg[i]);
    //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 %3.2f+-%3.2f\n",
	   Eg[i],Nphi_3g[i],eff_phi_3g[i],Nphi_3pi[i],eff_phi_3pi[i],
	   denom,xsec_3g[i],dxsec_3g[i],xsec_3pi[i],dxsec_3pi[i],
	   xsec_kpkm[i],dxsec_kpkm[i]);
  }

  TH2F *h=new TH2F("h","#phi cross section estimates",45,2.8,11.8,10,0,1.0);
  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_3g,dEg,dxsec_3g);
  g1->SetMarkerColor(8);
  g1->SetLineColor(8);
  g1->SetMarkerStyle(20);
  g1->Draw("p"); 
  
  TGraphErrors *g2=new TGraphErrors(45,Eg,xsec_3pi,dEg,dxsec_3pi);
  g2->SetMarkerColor(2);
  g2->SetLineColor(2);
  g2->SetMarkerStyle(20);
  g2->Draw("p");
  
  
  TGraph *g3=new TGraph(45,Eg,xsec_param);
  g3->Draw("l"); 
  TGraph *g4=new TGraph(45,Eg,xsec_param_up);
  g4->SetLineStyle(4);
  //g4->Draw("l"); 
  TGraph *g5=new TGraph(45,Eg,xsec_param_down);
  g5->SetLineStyle(4);
  //g5->Draw("l"); 

  
  TGraphErrors *g6=new TGraphErrors(45,Eg,xsec_kpkm,dEg,dxsec_kpkm);
  g6->SetMarkerColor(4);
  g6->SetLineColor(4);
  g6->SetMarkerStyle(24);
  g6->Draw("p");
  
  double oldxsec[1]={0.55};
  double doldxsec[1]={0.07};
  double Eold[1]={9.3};
  double dEold[1]={0.6};

  TGraphErrors *g7=new TGraphErrors(1,Eold,oldxsec,dEold,doldxsec);
  g7->SetMarkerStyle(20);
  g7->Draw("p");


  TLegend *legend=new TLegend(0.35,0.3,0.85,0.15);
  legend->AddEntry(g1,"#phi#rightarrow#eta#gamma","p");
  legend->AddEntry(g2,"#phi#rightarrow#pi^{+}#pi^{-}#pi^{0}","p"); 
  legend->AddEntry(g6,"#phi#rightarrowK^{+}K^{-}","p");
  legend->AddEntry(g7,"Old measurement, #phi#rightarrowK^{+}K^{-} (Ballam,PRD7,3150)","p");
  legend->AddEntry(g3,"#sigma parametrization (old data)","l");
  legend->Draw();
   
}