double myfunc(double *x,double *par){
  double m=x[0];  
  double m_pi=0.13957;  
  double f=0;
  double thresh=m_pi+m_pi;
  if (m>thresh){
    f+=par[0]*pow(m-thresh,par[1])*exp(1.+par[2]*m+par[3]*m*m);

    double m_rho=par[5];
    // m_rho=0.76;
    double a=m*m-m_rho*m_rho;
    double q2=0.25*(m*m-4.*m_pi*m_pi);
    double q2_rho=0.25*(m_rho*m_rho-4.*m_pi*m_pi);
    double R2=par[6];
    //par[7]=0.15;
    double Gamma_rho=par[7]*pow(q2/q2_rho,1.5)*m_rho/m*(1.+R2*q2_rho)/(1.+R2*q2);
    //Gamma_rho*=0.8;
    double b=m_rho*Gamma_rho;
    double b2=b*b;
    double a2=a*a;
    f+=par[4]*b2/q2/(a2+b2);
  
    double m_omega=par[9];
    m_omega=0.783;
  //par[4]=2.;
    double R2o=par[10]; 
    double ao=m*m-m_omega*m_omega;
    double q2_omega=0.25*(m_omega*m_omega-4.*m_pi*m_pi);
    double Gamma_omega=par[11]*pow(q2/q2_omega,1.5)*m_omega/m*(1.+R2o*q2_omega)/(1.+R2o*q2);
    //  Gamma_omega*=4.;
    double bo=m_omega*Gamma_omega;
    double b2o=bo*bo;
    double a2o=ao*ao;
    //  par[8]=0.;
    f+=par[8]*b2o/q2/(a2o+b2o);
    f+=2.*sqrt(par[4]*par[8])*b*bo/q2/(a2+b2)/(a2o+b2o)*(cos(par[12])*(a*ao-b*bo)+sin(par[12])*(a*bo-ao*b)); 
    
    
    double m_f2=par[14];
    double af=m*m-m_f2*m_f2;
    double q2_f2=0.25*(m_f2*m_f2-4.*m_pi*m_pi);
    double z=q2*par[15];
    double z0=q2_f2*par[15];
    double Gamma_f2=par[16]*pow(q2/q2_f2,2.5)*m_f2/m
      *((z0-3)*(z0-3)+9*z0)/((z-3)*(z-3)+9*z);
    double bf=m_f2*Gamma_f2;
    double b2f=bf*bf;
    double a2f=af*af;
    f+=par[13]*b2f/q2/(a2f+b2f);

   
  }

  return f;
}

void calc_rho_yield(int mybin,double &yield,double &unc){
  TF1 *f1=new TF1("f1",myfunc,0.28,2.,17);
  f1->SetParLimits(0,0.,1e6);
  f1->SetParameters(1,3.,1.8,-2.);
  f1->SetParameter(4,2000.);
  f1->SetParameter(5,0.76);
  f1->SetParameter(6,25.);
  f1->SetParLimits(6,0.,1e6);
  f1->SetParameter(7,0.15);
  f1->SetParameter(8,0.1);
  f1->SetParameter(9,0.78);
  f1->SetParameter(10,25.);
  f1->SetParLimits(10,0.,1e6);
  f1->SetParameter(11,0.008);
  f1->SetParameter(12,1.);
  f1->SetParameter(13,15);
  f1->SetParameter(14,1.27);
  f1->SetParameter(15,1.4);
  f1->SetParLimits(15,0.,1e6);
  f1->SetParameter(16,0.2);

  f1->SetRange(0.28,2.5);

  TH2F *hpippim=(TH2F*)_file0->FindObjectAny("PipPimMass_vs_Ebeam");
  TCanvas *c1=new TCanvas("c1","c1");

  char name[80];
  sprintf(name,"slice%d",mybin);
    
  
  TH1D *slice=hpippim->ProjectionY(name,mybin,mybin); 
  slice->GetXaxis()->SetRangeUser(0.28,2.5);
  slice->Fit("f1","r");
  
  char myname[80];
  sprintf(myname,"#pi^{+}#pi^{-} yield for E_{#gamma}=%3.1f [GeV]",2.9+0.2*(mybin-1));
  slice->SetTitle(myname);
  slice->SetYTitle("counts / 0.0025 GeV");
  slice->SetMinimum(0.);
  slice->Draw();
  TF1 *f2(f1);
  f2->SetParameter(4,0);
  f2->SetLineStyle(2);
  f2->Draw("same");
  
  double yield_unc=0;
  yield=slice->IntegralAndError(113,480,yield_unc)-f2->Integral(0.28,1.2)/0.0025;
  //yield=slice->Integral(201,400);//-f2->Integral(0.28,1.2)/0.0025;
  unc=sqrt(yield_unc*yield_unc+f2->Integral(0.28,1.2)/0.0025);


  cout <<"Rho:   " << yield << "+/-" << unc << endl; 


}