void calc_eta_xsecs(){
  double Eg[24];
  double dEg[24];

  double eff_eta_3piq[24]; 
  double deff_eta_3piq[24];
  double eff_eta_2g[24];
  double deff_eta_2g[24];

  double Nps[24];
  double Neta_2g[24],dNeta_2g[24];
  double Neta_3piq[24],dNeta_3piq[24];

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

  for (int i=0;i<24;i++){
    Eg[i]=7.1+0.2*double(i);
    dEg[i]=0.1;
  }

  ifstream eff_2g_file("eta_2g_eff.dat");
  for (int i=0;i<24;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<24;i++){
    eff_3pi_file >> eff_eta_3piq[i]; 
    eff_3pi_file >> deff_eta_3piq[i];
  }
  eff_3pi_file.close();


  ifstream in("eta_2g_yields.dat");
  for (int i=0;i<24;i++){
    double 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<24;i++){
    in2 >> Neta_3piq[i];
    in2 >> dNeta_3piq[i];
  }
  in2.close(); 

  double xsec_2g[24],xsec_3piq[24];
  double xsec_param[24],xsec_param_up[24],xsec_param_down[24];
  double dxsec_2g[24],dxsec_3piq[24];

  for (int i=0;i<24;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);
    // 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);
    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;
    /*
    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",24,7.0,11.8,10,0,0.4);
  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(22,Eg,xsec_2g,dEg,dxsec_2g);

  g1->SetMarkerStyle(20);
  g1->Draw("p"); 
 
  TGraphErrors *g2=new TGraphErrors(22,Eg,xsec_3piq,dEg,dxsec_3piq);
  g2->SetMarkerColor(2);
  g2->SetLineColor(2);
  g2->SetMarkerStyle(20);
  g2->Draw("p");
  
  
  //TGraph *g3=new TGraph(24,Eg,xsec_param);
  //g3->Draw("l"); 
  /*
  TGraph *g4=new TGraph(24,Eg,xsec_param_up);
  g4->SetLineStyle(2);
  g4->Draw("l"); 
  TGraph *g5=new TGraph(24,Eg,xsec_param_down);
  g5->SetLineStyle(2);
  g5->Draw("l"); 
  */

  TLegend *legend=new TLegend(0.45,0.7,0.85,0.85);
  legend->AddEntry(g1,"#eta#rightarrow2#gamma","p");
  legend->AddEntry(g2,"#eta#rightarrow#pi^{+}#pi^{-}#pi^{0}","p");
  //legend->AddEntry(g3,"#sigma parametrization (old data)","l");
  legend->Draw();

}