double my6gfunc(double *x,double *par){
  double f=0.;  

  //eta
  f+=par[0]*TMath::Gaus(x[0],par[1],par[2]);
  //eta-prime
  f+=par[3]*TMath::Gaus(x[0],par[4],par[5]);
  //f1
  f+=par[6]*TMath::Voigt(x[0]-par[7],par[2],par[8]);
  // background
  f+=par[9]*TMath::Gaus(x[0],par[10],par[11]);
  return f;
}

void calc_6g_yield(int mybin,double yields[],double uncs[]){
  TH2F *h6g=(TH2F*)_file0->FindObjectAny("SixGammaMass_vs_Ebeam");

  TF1 *f1=new TF1("f1",my6gfunc,0.3,2.5,12);
 
  TCanvas *c1=new TCanvas("c1","c1");

  char name[80];
  sprintf(name,"slice%d",mybin);
    
  
  TH1D *slice=h6g->ProjectionY(name,mybin,mybin);
  slice->GetXaxis()->SetRangeUser(0.5,0.6);
  slice->Fit("gaus","q");  
  if (slice->GetFunction("gaus")){
    f1->SetParameter(0,slice->GetFunction("gaus")->GetParameter(0));
    f1->SetParameter(1,0.55);
    f1->SetParameter(2,slice->GetFunction("gaus")->GetParameter(2)); 
  }

  slice->GetXaxis()->SetRangeUser(0.91,1.01);
  slice->Fit("gaus","q"); 
  f1->SetParameter(4,0.96);
  if (slice->GetFunction("gaus")){
    f1->SetParameter(3,slice->GetFunction("gaus")->GetParameter(0));
    f1->SetParameter(5,slice->GetFunction("gaus")->GetParameter(2));
  }

  slice->GetXaxis()->SetRangeUser(1.15,2.5);
  slice->Fit("gaus","q");
  if (slice->GetFunction("gaus")){
    f1->SetParameter(9,slice->GetFunction("gaus")->GetParameter(0));
    f1->SetParameter(10,slice->GetFunction("gaus")->GetParameter(1)); 
    f1->SetParameter(11,slice->GetFunction("gaus")->GetParameter(1));
  }
  f1->SetParameter(6,0.1*slice->GetBinContent(544));
  f1->SetParameter(7,1.3);
  f1->SetParameter(8,0.05);

  slice->GetXaxis()->SetRangeUser(0.3,1.6);
  //f1->Draw();
  slice->Fit("f1");
  
  TF1 *f2=new TF1("f2",my6gfunc,0.3,2.5,12);
  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));
  
  char myname[80];
  sprintf(myname,"6#gamma yield for E_{#gamma}=%3.1f [GeV]",2.9+0.2*(mybin-1));
  slice->SetTitle(myname);
  slice->SetYTitle("counts / 0.005 GeV");
  slice->SetMinimum(0.);
  slice->Draw();
  f2->SetLineColor(4);
  f2->Draw("same");
  //f2->Draw();

  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[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; 
  
 }