void calc_phi_t_xsecs2(double Eg){
  gStyle->SetOptStat(0);
  double eff[50],deff[50],eff1[50],deff1[50];
  double y[50],dy[50],y1[50],dy1[50];
  double t[50],dt[50];
  double Nps[42];
  double Ebincenter=Eg-0.1;
  double costheta[2040],xclas[2040],dxclas[2040],w[2040];
  double tclas[24];
  double xtclas[24];
  double dxtclas[24];
  
  double m_phi=1.01946;
  double m_phi_sq=m_phi*m_phi;
  double m_p=0.93827;
  double m_p_sq=m_p*m_p;
  unsigned int ind=0;
  ifstream inclas("clasdb_E63M1.txt");
  for (int i=0;i<2040;i++){
    inclas >> w[i];
    inclas >> costheta[i];
    inclas >> xclas[i];
    inclas >> dxclas[i];
    
    double s=w[i]*w[i];
    double E_phi=(s-m_p_sq+m_phi_sq)/(2*w[i]);
    double p_phi_sq=E_phi*E_phi-m_phi_sq;
    double p_phi=sqrt(p_phi_sq);
    double E_phi_sq=E_phi*E_phi;
    double p_beam=(s-m_p_sq)/2./w[i];
    double Eclas=(s-m_p_sq)/(2.*m_p);  
    if (Eclas>Eg+0.01) break;
    double myt=pow(E_phi-p_beam,2)-pow(p_phi-p_beam,2)
      -2.*p_phi*p_beam*(1.-costheta[i]);
    if (Eclas>Eg){
      double scale=0.5/(p_beam*p_phi);
      tclas[ind]=-myt;
      xtclas[ind]=scale*xclas[i];
      dxtclas[ind]=scale*dxclas[i];
      printf("%f %f\n",0.5/(p_beam*p_phi),2*s/(s-m_p_sq)/sqrt(s*s-2.*s*(m_phi_sq+m_p_sq)+pow(m_phi_sq-m_p_sq,2)));
      printf("%d E %f t %f dsigma/dt %f\n",ind,(s-m_p_sq)/(2.*m_p),myt,0.5*xclas[i]/(p_beam*sqrt(p_phi_sq)));
      ind++;
    }

  }


  int bin=int((Ebincenter-2.8)/0.2);
  int bin1=bin+1;

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

  double Nps_sum=Nps[bin];
  Nps_sum+=Nps[bin1];
  
  ifstream in1("phi_t_yields.dat");
  for (int i=0;i<bin;i++){
    string skipline;
    getline(in1,skipline);
    //   cout << skipline << endl;    
  }
  double dummy;
  in1 >> dummy;
  for (int i=0;i<50;i++){
    in1 >> y[i];
    in1 >> dy[i];
  }
  in1 >> dummy;
  for (int i=0;i<50;i++){
    in1 >> y1[i];
    in1 >> dy1[i];
  }
  
  in1.close();
  ifstream in2("phi_kk_t_acc.dat"); 
  for (int i=0;i<bin;i++){
    string skipline;
    getline(in2,skipline);
    //   cout << skipline << endl;    
  }
  in2 >> dummy;
  for (int i=0;i<50;i++){
    in2 >> eff[i];
    in2 >> deff[i];
  }  
  in2 >> dummy;
  for (int i=0;i<50;i++){
    in2 >> eff1[i];
    in2 >> deff1[i];
  } 

  in2.close();

  double y_acc[50],dy_acc[50];
  for (int i=0;i<50;i++){
    cout << eff[i] << endl;
    if (eff[i]>0.&&eff1[i]>0){
      double scale=1e6/1.293/Nps_sum/0.492/0.04;
      y_acc[i]=scale*(y[i]/eff[i]+y1[i]/eff1[i]);
      cout << y[i] << " " << y_acc[i] << endl;
      double eff_sq=eff[i]*eff[i];
      double eff1_sq=eff1[i]*eff1[i];
      dy_acc[i]
	=scale*sqrt(1./eff_sq*(dy[i]*dy[i]+y[i]*y[i]*deff[i]*deff[i]/eff_sq)
		    +1./eff1_sq*(dy1[i]*dy1[i]+y1[i]*y1[i]*deff1[i]*deff1[i]/eff1_sq));
			       
    }
    t[i]=0.02+0.04*i;
    dt[i]=0.02;
  }

  TCanvas *c2=new TCanvas("c2","c2");
  c2->Draw();
  TH2F *h=new TH2F("h","Differential cross section for #phi",10,0,1.5,10,0,2);
  h->SetXTitle("-t [GeV^{2}]");
  h->SetYTitle("d#sigma/dt [#mub/GeV^{2}]");
  h->Draw();

  TGraphErrors *g1=new TGraphErrors(46,&t[2],&y_acc[2],&dt[2],&dy_acc[2]);
  g1->SetLineColor(2);
  g1->SetMarkerColor(2);
  g1->SetMarkerStyle(20);
  g1->Fit("expo");
  g1->Draw("p");

  double ballamx[5]={1.5,0.68,0.23,0.15,0.04};
  double ballamdx[5]={0.3,0.15,0.11,0.10,0.03};
  double ballamt[5]={0.11,0.3,0.5,0.7,1.2};
  double ballamdt[5]={0.09,0.1,0.1,0.1,0.4};
 
  double ballamx_lo[5]={1.22,0.44,0.17,0.24,0.03};
  double ballamdx_lo[5]={0.22,0.12,0.10,0.10,0.02};
  double ballamt_lo[5]={0.1225,0.3,0.5,0.7,1.2};
  double ballamdt_lo[5]={0.0775,0.1,0.1,0.1,0.4};

  TGraphErrors *g2=NULL;
  if (Eg>5){
    g2=new TGraphErrors(5,ballamt,ballamx,ballamdt,ballamdx);
  }
  else{
    g2=new TGraphErrors(5,ballamt_lo,ballamx_lo,ballamdt_lo,ballamdx_lo);
  }
  g2->SetLineColor(8);
  g2->SetMarkerColor(8);
  g2->SetMarkerStyle(22);
  g2->Draw("p");

  TGraphErrors *g3=new TGraphErrors(24,tclas,xtclas,0,dxtclas);
  g3->SetLineColor(4);
  g3->SetMarkerColor(4);
  g3->SetMarkerStyle(25);
  g3->Draw("p");

  TLegend *legend =new TLegend(0.4,0.5,0.8,0.7);
  char mytitle[80];
  sprintf(mytitle,"our measurement at E#gamma=%3.1f GeV",Eg);
  legend->AddEntry(g1,mytitle,"p");
  legend->AddEntry(g2,"Ballam 1973","p"); 
  legend->AddEntry(g3,"CLAS 2013","p");
  legend->Draw();


}