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

  double eff_rho[45]; 
  double deff_rho[45];
  double Nps[45];
  double Nrho[45],dNrho[45];

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

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

  ifstream eff_2pi_file("rho_eff.dat");
  for (int i=0;i<45;i++){
    double dummy;
    eff_2pi_file >> dummy;
    eff_2pi_file >> eff_rho[i]; 
    eff_2pi_file >> deff_rho[i];
  }
  eff_2pi_file.close();


  ifstream in("rho_yields.dat");
  for (int i=0;i<45;i++){
    in >> Eg[i];
    in >> Nrho[i];
    in >> dNrho[i];
  }
  in.close(); 

  double xsec_rho[45];
  double xsec_param[45],xsec_param_up[45],xsec_param_down[45];
  double dxsec_rho[45];

  for (int i=0;i<45;i++){
    double denom=Nps[i];
    double rho_numer=Nrho[i];
    double r_rho=rho_numer/denom;
    double dr_rho=0.;
    if (rho_numer>0)
      dr_rho=r_rho*sqrt(1./denom+(dNrho[i]*dNrho[i])
		       /rho_numer/rho_numer);
    double scale_rho=1e6/(1.293)/eff_rho[i];

    xsec_rho[i]=scale_rho*r_rho;
    dxsec_rho[i]=scale_rho*dr_rho;

    xsec_param[i]=1000.*(0.813e-2+0.0192*pow(Eg[i],-0.59));
    double xsec_var=0.357e-3*0.357e-3+0.113e-2*0.113e-2*pow(Eg[i],-2.*0.59)
      +0.192e-1*0.192e-1*pow(Eg[i],2.*(-1.59))*0.27*0.27*0.59*0.59;
    //printf("%f\n",sqrt(xsec_var));
    xsec_param_up[i]=xsec_param[i]+1000.*sqrt(xsec_var);
    xsec_param_down[i]=xsec_param[i]-1000.*sqrt(xsec_var);

    printf("%3.2f %f %f %f %3.2f+-%3.2f\n",
	   Eg[i],Nrho[i],eff_rho[i],
	   denom,xsec_rho[i],dxsec_rho[i] );
  }

  TH2F *h=new TH2F("h","#rho cross section estimates",45,2.8,11.8,10,0,30.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 *g2=new TGraphErrors(45,Eg,xsec_rho,dEg,dxsec_rho);
  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(g2,"#rho#rightarrow#pi^{+}#pi^{-}","p");
  legend->AddEntry(g3,"#sigma parametrization (old data)","l");
  legend->Draw();

}