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

  double eff_omega_3pi[45]; 
  double deff_omega_3pi[45];
  double eff_omega_3g[45];
  double deff_omega_3g[45];

  double Nps[45];
  double Nomega_3g[45],dNomega_3g[45];
  double Nomega_3pi[45],dNomega_3pi[45];

  //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_3g_file("omega_3g_eff.dat");
  for (int i=0;i<45;i++){
    eff_3g_file >> eff_omega_3g[i];
    eff_3g_file >> deff_omega_3g[i];
  }
  eff_3g_file.close(); 

  ifstream eff_3pi_file("omega_3pi_eff.dat");
  for (int i=0;i<45;i++){
    eff_3pi_file >> eff_omega_3pi[i]; 
    eff_3pi_file >> deff_omega_3pi[i];
  }
  eff_3pi_file.close();


  ifstream in("omega_yields.dat");
  for (int i=0;i<45;i++){
    double dummy;
    in >> dummy;
    in >> Nomega_3g[i];
    in >> dNomega_3g[i];
    in >> Nomega_3pi[i];
    in >> dNomega_3pi[i];
  }
  in.close(); 

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

  for (int i=0;i<45;i++){
    double denom=Nps[i];
    double omega_3g_numer=Nomega_3g[i];
    double r_omega_3g=omega_3g_numer/denom;
    double dr_omega_3g=0.;
    if (omega_3g_numer>0)
      dr_omega_3g=r_omega_3g*sqrt(1./denom+(dNomega_3g[i]*dNomega_3g[i])
		       /omega_3g_numer/omega_3g_numer);
    // printf("%f %f %f \n",1/denom,dNomega_3g[i],dNomega_3g_acc[i]);
    
    double scale_omega_3g=1e6/(1.293)/eff_omega_3g[i]
      /0.084;
			     //myflux << denom*1.249/(eff_ps[i]*1.65e-4) << endl;
    
    xsec_3g[i]=scale_omega_3g*r_omega_3g;
    dxsec_3g[i]=scale_omega_3g*dr_omega_3g;

    double omega_3pi_numer=Nomega_3pi[i];
    double r_omega_3pi=omega_3pi_numer/denom;
    double dr_omega_3pi=0.;
    if (omega_3pi_numer>0)
      dr_omega_3pi=r_omega_3pi*sqrt(1./denom+(dNomega_3pi[i]*dNomega_3pi[i])
		       /omega_3pi_numer/omega_3pi_numer);
    double scale_omega_3pi=1e6/(1.293)/eff_omega_3pi[i]
      /0.892;

    xsec_3pi[i]=scale_omega_3pi*r_omega_3pi;
    dxsec_3pi[i]=scale_omega_3pi*dr_omega_3pi;

    xsec_param[i]=1.14+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],Nomega_3g[i],eff_omega_3g[i],Nomega_3pi[i],eff_omega_3pi[i],
	   denom,xsec_3g[i],dxsec_3g[i],xsec_3pi[i],dxsec_3pi[i] );
  }

  TH2F *h=new TH2F("h","#omega cross section estimates",45,2.8,11.8,10,0,4.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->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(2);
  g4->Draw("l"); 
  TGraph *g5=new TGraph(45,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,"#omega#rightarrow#pi^{0}#gamma","p");
  legend->AddEntry(g2,"#omega#rightarrow#pi^{+}#pi^{-}#pi^{0}","p");
  legend->AddEntry(g3,"#sigma parametrization (old data)","l");
  legend->Draw();

}