double my2gfunc(double *x,double *par){
  double f=0.;
  //eta
  f+=par[0]*TMath::Gaus(x[0],par[1],par[2]);
  f+=par[3]*TMath::Gaus(x[0],par[4],par[5]);
  //omega
  f+=par[6]*TMath::Gaus(x[0],par[7],par[8]);
  f+=par[9]*TMath::Gaus(x[0],par[10],par[11]);
  //eta-prime
  f+=par[12]*TMath::Gaus(x[0],par[13],par[14]);
  // higher mass
  f+=par[15]*TMath::Gaus(x[0],par[16],par[17]);
  // background
  f+=par[18]+par[19]*x[0]+par[20]*x[0]*x[0];
  return f;
}

void calc_2g_yield(int mybin,double yields[],double uncs[]){
  TH2F *h2g=(TH2F*)_file0->FindObjectAny("TwoGammaMass_vs_Ebeam");

  TF1 *f1=new TF1("f1",my2gfunc,0.3,2.,21);
 
  // TCanvas *c1=new TCanvas("c1","c1");

  char name[80];
  sprintf(name,"slice%d",mybin);
  //  double yields[3],uncs[3];
    
  
  TH1D *slice=h2g->ProjectionY(name,mybin,mybin);
  slice->GetXaxis()->SetRangeUser(0.525,0.575);
  slice->Fit("gaus","q");
  if (slice->GetFunction("gaus")){
    double height=slice->GetFunction("gaus")->GetParameter(0);
    f1->SetParameter(0,height);
    f1->SetParameter(3,0.1*height);
  }
  f1->SetParameter(1,0.55);
  f1->SetParameter(2,0.01);
  f1->SetParameter(4,0.55);
  f1->SetParameter(5,0.02); 
  slice->GetXaxis()->SetRangeUser(0.7,0.85);
  slice->Fit("gaus","q"); 
  if (slice->GetFunction("gaus")){
    double height=slice->GetFunction("gaus")->GetParameter(0);
    f1->SetParameter(6,height); 
    f1->SetParameter(7,0.77);
    f1->SetParameter(8,0.03);
    f1->SetParameter(9,0.5*height);
    f1->SetParameter(10,0.73);
    f1->SetParameter(11,0.03);
  } 
  slice->GetXaxis()->SetRangeUser(0.92,1.);
  slice->Fit("gaus","q"); 
  if (slice->GetFunction("gaus")){
    double height=slice->GetFunction("gaus")->GetParameter(0);
    f1->SetParameter(12,height); 
    f1->SetParameter(13,0.96);
    f1->SetParameter(14,0.01);
  } 
  slice->GetXaxis()->SetRangeUser(1.1,1.5);
  slice->Fit("gaus","q"); 
  if (slice->GetFunction("gaus")){
    double height=slice->GetFunction("gaus")->GetParameter(0);
    f1->SetParameter(15,height); 
    f1->SetParameter(16,slice->GetFunction("gaus")->GetParameter(1));
    f1->SetParameter(17,slice->GetFunction("gaus")->GetParameter(2));
  }

  slice->GetXaxis()->SetRangeUser(0.3,2.);
  slice->Fit("f1");
  
  TF1 *f2=new TF1("f2",my2gfunc,0.3,2.,21);
  f2->SetParameter(6,f1->GetParameter(6));
  f2->SetParameter(7,f1->GetParameter(7));
  f2->SetParameter(8,f1->GetParameter(8)); 
  f2->SetParameter(9,f1->GetParameter(9));
  f2->SetParameter(10,f1->GetParameter(10));
  f2->SetParameter(11,f1->GetParameter(11)); 
  f2->SetParameter(15,f1->GetParameter(15));
  f2->SetParameter(16,f1->GetParameter(16));
  f2->SetParameter(17,f1->GetParameter(17)); 
  f2->SetParameter(18,f1->GetParameter(18));
  f2->SetParameter(19,f1->GetParameter(19));
  f2->SetParameter(20,f1->GetParameter(20));

  slice->GetXaxis()->SetRangeUser(0.3,2.);
  
  char myname[80];
  sprintf(myname,"2#gamma 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();
  f2->SetLineColor(4);
  f2->Draw("same");

  double yield_unc=0.;
  yields[0]=slice->IntegralAndError(201,240,yield_unc)-f2->Integral(0.5,0.6)/0.0025; 
  uncs[0]=sqrt(yield_unc*yield_unc+f2->Integral(0.5,0.6)/0.0025); 
  // yields[2]=f4->Integral(0.96,1.08)/0.0025; 
  yields[1]=slice->IntegralAndError(365,404,yield_unc)-f2->Integral(0.91,1.01)/0.0025;
  uncs[1]=sqrt(yield_unc*yield_unc+f2->Integral(0.91,1.01)/0.0025);
  cout <<"Eta: " << yields[0]<< "+/-" << uncs[0]<< endl; 
  cout <<"Eta-prime:   " << yields[1]<< "+/-" << uncs[1] << endl; 

 }